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ABSTRACT 


An  equation  of  state  suitable  for  calculating  the 
compression  of  a  melting  solid  is  described.  Some  elementary 
ideas  about  melting  are  reviewed  and  some  standard  relations 
between  P  and  T  in  the  melting  region  are  described.  The  equa¬ 
tion  of  state  and  melting  law  are  combined  in  a  program  for 
calculating  the  Hugoniot  through  the  mixed  phase  region.  Results 
are  described  for  lead,  which  melts  at  a  shock  pressure  of  about 
400  kilobars  with  a  Kennedy  equation  and  700  kilobars  for  a 
Simon  equation. 

The  Eyring  theory  for  equation  of  state  of  liquids  is 
examined  for  argon,  and  Hugoniot  curves  are  calculated.  Cal¬ 
culations  agree  with  the  most  dense  case  of  van  Thiel  and  Alder 
to  13  kilobars,  then  depart  dramatically  from  measured  values. 

The  theory  of  plastic  wave  propagation  in  two -dimens ions 
is  discussed  and  calculations  of  allowed  directions  are  described. 
These  will  ultimately  be  of  use  in  discussing  the  reflection  of 
obliquely  incident  waves  in  an  elastic-plastic  medium. 

Some  of  the  basic  physical  mechanisms  in  solid- solid 
phase  transitions  are  reviewed  and  the  applicability  of  thermo¬ 
dynamics  to  such  transitions  is  brought  into  question.  An 
elementary  model  for  a  non-equilibrium  transition  in  iron  is 
suggested  and  p-v  calculations  are  made  for  several  values  of 
the  parameters.  It  is  evident  that  no  conclusions  about  the 
time  dependence  of  the  ot-t  transition  can  be  drawn  from  second 
state  shock  measurements,  although  it  may  be  possible  to  infer 
useful  information  about  metastable  states. 
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INTRODUCTION 


This  project,  since  its  inception,  has  been  concerned 
with  problems  of  equations  of  state,  constitutive  relations, 
phase  transitions,  and  wave  propagation.  These  are  large  prob¬ 
lems,  not  to  be  wrapped  up  and  put  aside  in  one-year  packages  by 
the  part-time  efforts  of  one  man,  or  of  several  men.  Instead, 
once  having  got  started  they  stick  in  the  mind  and  bits  and 
pieces  of  new  understanding  or  new  accomplishments  come  along, 
sometimes  in  unexpected  directions.  When  it  comes  time  to  put 
these  together  in  a  summary  report  at  the  end  of  a  year,  their 
unity  is  not  always  apparent;  it  is  not  always  clear  that  they 
are  parts  of  a  whole  which  isn't  very  easy  to  subdivide.  So  it 
is  with  the  present  report;  so  it  is  presented  in  three  distinct 
parts  as  it  was  worked  out  by  three  different  men  working  on  the 
problems  and  their  assistants. 

Part  A  deals  quite  directly  and  clearly  with  the  stated 
objectives  of  the  contract;  it  comprises  a  relatively  straight¬ 
forward  and  tedious  calculation  of  the  equilibrium  shock  Hugoniot 
of  a  melting  solid.  In  the  process  of  doing  this  calculation  the 
author  has  experienced  some  revealing  insights  into  the  features 
of  a  total  equation  of  state,  and  a  side  excursion  into  the  theory 
of  equations  of  state  of  liquids  has  shown  that  the  Eyring  Sig¬ 
nificant  Structure  Theory  may  be  amenable  to  modifications  which 
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would  lead  to  a  total  equation  of  state,  when  coupled  with 
suitable  data. 

Part  B  is  a  beginning,  an  introduction  to  the  quite  com¬ 
plex  subject  of  propagation  of  multiple  waves  in  anelastic, 
yielding  materials.  It  shows  even  at  this  beginning  stage  that 
the  mechanics  is  more  complicated  than  we  believed  when  we 
thought  only  of  plane  waves  following  parallel  plane  waves.  In 
its  present  stage  it  begins  to  provide  a  foundation  for  under¬ 
standing  wave  structures  in  more  realistic  solid  models  than  have 
been  heretofore  commonly  used.  It  may  even  eventually  provide 
for  better  interpretation  of  experiments  and  designs  for  experi¬ 
ments  . 

Part  G  is  something  quite  different — a  summarizing  of 
physical  ideas  about  the  causes  of  phase  transitions  and  the  mech 
ar.ics  of  their  occurrence.  The  metallurgist  has  long  been  aware 
that  equilibrium  thermodynamics  plays  only  a  minor  role  in  solid- 
solid  phase  transitions.  The  thermodynamicist  and  the  physicist 
are  just  beginning  to  learn  this.  This  summary  suggests  an 
elementary  model  for  martensitic  phase  transition.  A  calculation 
for  iron  shows  that  the  quasi-stable  p-v  relation  in  the  mixed 
phase  can  be  varied  almost  at  will  by  varying  the  assumed  metal¬ 
lurgical  parameters. 


PART  A 


SHOCK  PROPAGATION  AND  MELTING 
G.  E.  Duvall 

I.  Melting  Phase  Boundaries  in  the  P-V  Plane 

In  Fig.  1,  let  OCFGBJ  be  the  coexistence  region  for 
solid  and  liquid;  9-^  is  the  solid  phase  and  92  is  liquid. 
Suppose  ABCD  to  be  the  Hugoniot.  We  wish  to  determine  first  the 
phase  boundaries  GBJ  and  FCO  in  terms  of  known  quantities. 

The  transition  is  first  order,  and  at  equilibrium  the 
Clausius-Clapeyron  equation  obtains: 

dP/'dT  =  AS/AV  =  AH/TaV  »  f(P)  (1.1) 

AV  =  V2(P,T)  -  V1(P,T)  (1.2) 

The  rate  of  change  of  entropy  with  temperature  on  the  phase 
boundary  JBG  is 

d3x/dT  =  (Cpl/T)  -  (SV/dT)pl(dP/dT) 

-  Cpl(dP/dT)[(dT/TdP)  -  U-/Cpl)(SV/BT)plJ  (1.3) 

where  subscript  "pi"  denotes  a  quantity  evaluated  at  constant 
pressure,  in  phase  1,  on  the  phase  boundary  GBJ. 

We  assume  that  Cp-^  >  0  and  (&V/S T)p^  >  0  .  Then  if 
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If  d?/dT  >  0,  the  sign  of  dS^/dT  depends  upon  the  magnitude 
of  dP/dT;- 

dT/dP  >  T(av/3T)pl/Cpl 
dS1/dT  >  0 

dT/dP  <  T(9V/9T)pl/Cpl 
dS1/dT  <  0 

We  may  summarize  these  relations  as  in  Table  I.  defining  three 
distinct  types  of  phase  transitions. 


Table  I 

Classification  of  Phase  Transitions 


AV 

-  V,  -  vs  - 

V2  -  Vj  >  0 

Type  of 
Transition 

d3^/dT 

dP/dT 

AS 

I 

>  0 

<  0 

<  0 

II 

<  0 

>  0 

>  0 

III 

>  0 

>  0 

>  0 

In  terms  of  observables  at  the  melting  point,  the  defining 
conditions  are 

dP/dT  <  0  Type  I 

dP/dT  >  0;  AVm/AHm  <  (dV/3T)pl/Cpl  TyPe  11 
dP/dT  >  0;  AV^/AH^  >  (&V/3T)pl/C 


Type  III 
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where  Av^  =  -  Vgo^  =  volume  expansion  on  melting 

AHp  =  latent  heat  of  melting. 

In  this  way  melting  transitions  for  which  Av  >  0  are 
divided  into  three  exhaustive  categories.  The  Metals  Reference 
Handbook  shows  that  Type  III  prevails  among  the  metals  (Table  II) . 

Of  these  three  categories,  II  and  III  may  be  considered  normal, 
i.e.,  S  increases  on  melting. 

The  relative  slopes  of  isotherm,  adiabat,  and  phase 
boundary  may  be  found  as  follows: 

In  the  solid  phase, 

(av/ap)sl  =  (*v/ap)T1  +  (&v/aT)pl  (*T/ap)sl  (1.4) 

Also,  the  slope  of  the  phase  boundary  may  be  written  as 

dV1/dP  =  (av/ap)n  +  (*V/dT)pl  (dT/dP)  (1.5) 

For  transitions  of  Type  III: 

(9T/ap)sl  =  T(dV/BT)pl/Cpl  <  dT/dP  (1.6) 

if  ar*d  (AV/oT)p  are  both  positive. 

Substitution  of  this  inequality  into  Eq.  (1.5)  and  comparing  with 
Eq.  (1.4)  yields  the  result 

dVp/dP  >  (3V/3P)S1  (1.7) 

Adiabats  in  the  solid  phase  exit  from  the  coexistence  region 
with  negative  slope  and  increasing  pressure. 

A  closely  related  treatment  has  recently  appeared  in  the 
literature.  See  Ref.  1.  However,  the  author  assigns  lead  to  Type  II 
instead  of  Type  III  as  a  result  of  numerical  error. 


The  three  types  of  transition  can  also  be  characterized 
in  terms  of  the  relative  slopes  of  isotherms,  adiabats  and  the 
cylinder  defining  the  mixed  phase  region  in  P-V-T  space.  First 
construct  the  P-V-T  surface  for  the  solid,  extending  it  to 
arbitrary  T  in  a  metastable  state.  Suppose  the  mixed  phase 
cylinder  to  be  degenerate,  i.e.,  a  plane,  in  order  to  simplify 
the  description.  The  plane  always  lies  parallel  to  the  V-axis. 
When  it  is  also  perpendicular  to  the  T-axis,  its  intersection 
with  the  P-V-T  surface  of  the  solid,  i.e.  the  phase  boundary, 
coincides  with  an  isotherm  of  the  solid.  This  is  the  case 
dT/dP  =  0. 

If  the  plane  inclines  toward  smaller  T  as  P  increases, 
dT/dP  <  0  and  the  transition  is  of  Type  I.  If  it  inclines 
toward  larger  T,  it  describes  first  Type  II  in  which  the  phase 
boundary  splits  the  isotherms  and  adiabats  and  then  Type  III  in 
which  both  isotherms  and  adiabats  issue  from  the  phase  boundary 
with  increasing  P.  These  three  cases  are  illustrated  in 


T  (b) 

Fig. 3 

Phase  Boundary  for  Type  m  Melting 

(a)  P-V  Plane 

(b)  S-T  Plane 
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II.  Adiabats  in  the  Mixed  Phase  Region 

Adiabats  in  the  mixed  phase  region  are  shown  in  the 
P-V  plane  in  Fig.  4.  The  entropy  at  point  C  is: 

,,B 

SC  =  Sa  +  J  (dS1/dP)d?  +  x&S  (2.1) 

A 

where 

x  =  (V-V-,  )/(V0-V, )  =  fraction  of  mass  in  phase  2 

at  point  C  (2.2) 

AS  =  S(B')  -  S(B)  (2.3) 

V  =  V(C),  V2  =  V(B'),  Vx  =  V(B) 

Differentiating  Eq.  (2.1)  yields 

(dS1/dP)  +  (AS/AV)((dV/SP)SM  -  (dVx/dP)) 

+  (v-v^  d(AS/AV)/dP  =  0  (2.4) 

where  (dV/dP)SM  =  ac^ia^atic  derivative  in  the  mixed  phase 
region,  R.  Neglecting  the  last  term  in  Eq.  (2.4)  yields  an 
expression  for  the  adiabat 

(&v/d?)SM  -  (dVL/dP)  -  (dS1/dP)(dT/dP) 

=  (BV/SP)T1  +  2  (av/3T)pl(dT/dP)  -  Cpl(dT/dP) 2/T  (2.5) 

The  discontinuity  in  the  slope  of  the  adiabat  is 

(3V/SP)S1  -  <av/ap>SM1  =  [(T/Cpl)%  (8V]/8T)p  -  (cpl/T)% 

(dT/dP)J2  >  0  ; 


(2.6) 


DISCONTINUITY  IN  ADIABAT  AT  PHASE  BOUNDARY  FOR  TYPE  3 E 

TRANSITION 
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i.e.,  the  adiabat  is  always  discontinuous  at  the  phase  boundary 
and  the  sign  of  the  discontinuity  is  such  that 

|  ap/av  lsl>  |  ap/av  |SM1  .  (2.7) 

III.  Construction  of  the  Hugoniot 

We  suppose  the  material  to  be  shocked  from  a  point  in 
the  solid  state  to  a  final  state  which  may  lie  in  the  solid,  in 
the  mixed  phase  region,  or  in  the  liquid  state.  It  is  not 
known  in  advance  whether  the  final  state  is  reached  through 
single  or  multiple  shocks,  consequently  it  is  most  appropriate 
to  construct  the  Hugoniot  incrementally,  examining  at  each  step 
to  determine  whether  or  not  a  new  shock  is  initiated. 

A.  Solid  Phase  Hugoniot 

The  differential  equation  of  the  Hugoniot  curve  for  an 
equation  of  state  of  the  form  of  Eq.  (4.1)  in  a  single  phase  is  (1) 

dP/dv  =  [<ap/av)s  +  <r/2v)(p  -  pq)]/[i  -  (r/2v)(vo2  -  v)j  (3.1) 

where  (PQ,  V  is  the  initial  state  and  r  is  the  Gruneisen 
ratio.  We  now  consider  whether  or  not  a  single,  stable  shock 
to  pressure  P^  will  also  be  stable  to  P^  +  &P  .  Since  a 
single  shock  from  PQ  to  P  is  assumed  to  be  stable,  the  Rayleigh 
line  connecting  (P  ,  V  )  with  (Pp  V^)  intersects  the  Hugoniot 
curve  only  at  those  two  points; 

Numbers  in  parentheses  reference  to  references  on 
pp.  39-40. 


-  dP/dV) 


2  (Pl"Po)/(Vo2“Vl) 


(3.2) 


RH,P1 


Here  P^  denotes  the  slope  on  the  lower  side  of  the  point 
Pp  if  the  Hugoniot  is  discontinuous  there.  If  the  point 
P1  +  6P1  a-'-so  t0  be  attained  through  a  single  shock, 
Condition  (3.2)  must  hold  on  the  upper  side  of  P^,  denoted 
by  Px+  : 


-  dP/dVW  *  (PrPo)/(Vo2'Vl)  <3-3> 

Substitution  of  Eq.  (3.1)  into  (3.3)  yields  the  condition  for 
stability: 


lajVtJj-u/  -  r^-vp/zv^  /  L1  -  (ri/2V1)(Vo2-V1)J 


*  1  (3.4) 


where 

4  -  V  <sp/av>SjPi+ 

(urui)2  -  vi2  (prpo)/(vo2'vi)  • 

If  r(V0?Vi)/2Vi  <  1  the  stability  condition  further  reduces  to 

a? / (U-u-i ) 2)  ,  2  1  .  (3.5) 

1  1  PI 

B.  Mixed  Phase  Hugoniot 

In  Fig.  6,  the  region  ABCD  denotes  the  part  of  the  mixed 
phase  region  through  which  the  Hugoniot  passes.  F  is  the  inter¬ 
section  of  the  solid  phase  Hugoniot  with  the  boundary  between 


Hugoniot  in  the  Mixed  Phase  Reg 


solid  and  mixed  phase.  G  is  a  point  on  the  mixed-phase 
equilibrium  Hugoniot.  The  enthalpy  difference  between  G  and 
F  can  be  written 

rJ 

Hg  -  Hf  =  J  (dH/dP)d?  !  ^(V-V1)/(V2-V1)  (3. 

F 

where  dH/dP  is  the  variation  of  enthalpy  along  the  phase 
boundary  and  AH(P,T)  is  the  enthalpy  difference  between  liquid 
and  solid  at  pressure  P,  temperature  T: 

AU  =  TA7  -T7  \  JTD/J'T'  —  _  tt  /  O 

JL  \ » 2  vi x  /  ux  n.j^  •  \  J  • 

The  Clausius-Clapeyron  equation  has  been  used  to  obtain  this 
result.  Substitution  of  Eq.  (3.7)  into  (3.6)  yields  the  result 

Hg  -  Hp  =  Hj  -  Hy  +  T(V-V1)dP/dT  .  (3. 

From  the  Rankine-Hugoniot  relation  we  have,  for  a  single  shock 
from  0  to  G: 

HG  -  Ho2  =  (P-?0)(Vo2+V)/2  (3. 

Combining  Eqs.  (3.8)  and  (3.9)  yields  the  result: 

(P-P0)(Vc2+V)/2  =  +  T(V-V1)(dP/dT)  -  Hq2  (3. 

Differentiate  Eq.  (3.10)  to  obtain: 

dF/dV  =  ((TdP/dT)  -  (P-P0)/2)/(A'  +B')  (3. 

A'  =»  (Vq2“V)/2  -  T(V-V1)(dT/dP)(d2P/dT2)  (3. 


B'  = 


V,  -  dH, /dP  +  T  dV, /dT 


(3. 


V1  -  dH^/ dP  =  T(3V1/5T)p  -  Cpl  dT/dP 

where  subscript  "1"  denotes  quantities  evaluated  in  the  solid 


(3.14) 


phase  at  the  phase  boundary.  Two  further  relations  are 
required: 

TdV1/dT  =  T(3V1/ST)p  +  T(av1/5P)T(dP/dT)  (3.15) 

cPi  =  Si  +  T(5pi/ST)v  (avi/ST)p  (3.15) 

Substitution  of  Eqs.  (3.14)-(3.16)  into  (3 . 13)  yields : 

B'  =  2T(Bv1/ST)p  -  CV1  dT/dP  -  T(aP1/ST)y  (Bv1/ST)p  (dT/dP) 

(3.17) 

Combining  this  with  Eq.  (3.12)  yields 
A'  +  B'  =  B  -  A 

where 

B  =  2T(BV1/dT)p  +  T(BV^/dp)T  (dP/dT)  +  (Vo2-V)/2  (3.18) 

A  =  (cyl  +  T(dP1/dT)v  (<iV/dT)p  +  T(V-V1)(d2P/dT2)) (dT/dP) 

(3.19) 

Then 

dP/dV  =  ( (P-PQ) / 2  -  TdP/dT)/(A-B)  (3.20) 

Equation  (3.20)  is  equivalent  to  one  given  by  V.  D.  Urlin  and 
A.  A.  Ivanov  in  Ref.  2. 

C.  Liquid  Phase 

At  the  boundary  between  mixed  phase  and  liquid  phase 
a  test  must  be  made  according  to  Eq.  (3.3)  to  determine  whether 
or  not  a  single  shock  into  the  liquid  phase  is  stable.  If  it 
is,  the  integration  is  continued  to  higher  values  of  P,  using 
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the  generalized  form  of  Eq.  (3.1): 

dp/dv  =  (<ap/av)s  +  (p-Po)(ap/5T)v/2cv)  / 

(i  -  (Vo2-V)(SP/aT)v/2Cv)  .  (3.21) 

IV.  Equations  of  State 
A.  Solid  Phase 

Calculation  of  the  Hugoniot  curve  through  the  solid  and 
mixed  phases  can  be  accomplished  using  only  the  equation  of 
state  of  the  solid.  Since  both  the  Kugoniot  and  the  phase 
boundary  must  be  calculated  in  the  P~V  plane,  the  equation  of 
state  must  be  complete  and  internally  consistent.  For  simplicity 
we  choose  a  Mie-Gruneisen  equation  with  Debye  variations  of 
thermal  energy.  As  we  shall  see,  this  produces  some  disagreement 
with  measurements,  but  for  the  present  we  ignore  these  for  the 
sake  of  theoretical  consistency.  The  Mie-Gruneisen  equation  is 

?(V,E)  =  pk(v)  +  (r/v)  (e(v,t)  -  Ek(V))  (4.1) 

where  P^(V)  and  E^(V)  are  pressure  and  internal  energy 
on  the  0°K  isotherm,  and  F  is  the  Gruneisen  parameter: 

r(V)  =  (V/Cv)(9P/aT)v  (4.2) 

where  Cy  is  specific  heat  at  constant  volume.  Following 
Rice,  McQueen  and  Walsh  (13)  we  write: 

pK(V)  =  PK(lO  =  +  b^2  +  b3^3  +  PQ 


(4.3) 


>1 
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r(V)  =  aQ  +  a^n  +  a2n2  +  a3P-3 


V*  =  (VQ/V)  -  1  =  (P/PQ)  -  1 


V  =  specific  volume  at  P  =  P  ,  T  =  0°K 
o  o 

The  difference  E  -  E  Eth  is  the  thermal  energy: 

Eth  =  3RT  D(9/T) 

where  6  is  the  Debye  temperature,  R  is  the  gas  constant 
divided  by  the  molecular  weight,  and  D(6/T)  is  the  Debye 
function: 

/  on  >e/T  - 

D(9/T)  =  (3/(6/T)JJj  (xJ/(exp(x)  -  l))dx  . 

o 

For  small  values  of  9/T,  D  can  be  expanded  in  series: 


D 


1  -  .375  (0/T)  +  .05  (9/T)‘ 


This  is  accurate  to  .3%  for  6/t  <  .3 

The  Debye  temperature  is  related  to  the  Gruneisen 
parameter  by  the  equation 

d^nQ/d^nV  =  -T 

Using  Eq.  (4.4)  for  T,  this  integrates  to 


9  =  0O  exp((aQ-a1+a2-a3)  l n(i*+l) 


+  (a^-a2+a3)u  +  (a2~a3)  v-  /2 


+  a3P- 


V3) 


(4,4) 


(4.5) 


(4.6) 


(4.7) 


(4.8) 


(4.9) 
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The  specific  heat  is  defined  as 

Cy  =  (d  Eth/ST)V  (4.10) 

=  3r(4D  -  3(9/T)/(exp(8/T)  -  1))  (4.11) 

~  3r(i  -  .05(e/T)2)  .  (4.12) 

In  calculating  the  Hugoniot,  we  shall  need  some 
thermodynamic  coefficients  which  can  be  calculated  from  the 


above  equations: 

d(r/v)/dV  =  -  (VQ/V3) (a1+2a2P  +  3a3^2)  -  r/v2  (4.13) 

(9  Eth/3V)T  =  3R  D*  de/dV  (4.14) 

■  -  3Rre(-.375  +  .l(e/T))/V  (4.15) 

D'  =  dD/d(e/T)  (4.16) 

dPK/dV  =  (-VQ/V2)(b1  +  2b2^  +  3b3n2)  (4.17) 

(4P/dV)T  =  dPK/dV  +  Eth  d(T/V)/dV  +  (I*/V)  (dEth/*V)T  (4.18) 
(av/aT)p  =  -  rcv(sv/ap)T/v  (4.19) 

(ap/dv)s  -  (ap/av)T  -  r2  cvt/v2  (4. 20) 


B.  Liquid  Phase 

One  of  the  most  promising  theories  of  the  liquid  state 
for  computational  purposes  is  Henry  Eyring's  Significant  Structure 
Theory  (Ref.  9) .  The  essence  of  the  theory  is  that  a  liquid 
consists  of  a  solid  containing  holes  of  atomic  dimensions,  and 
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that  the  holes  behave  like  the  molecules  of  an  ideal  gas.  In 
order  to  test  it  at  high  presewps.  isotherms,  adiabats  and 
Hugoniot  curves  for  argon  have  been  calculated  from  it  and  the 
last  compared  with  measured  values  reported  in  the  literature 
(Ref.  12).  The  two  agree  remarkably  well  at  low  pressures, 
but  the  theory  fails  to  properly  account  for  the  energy  of  cold 
compression  and  the  computation  fails  at  higher  pressures  (>  12 
Kbars  for  the  highest  initial  density) .  These  computations  are 
reported  in  Appendix  B. 

Various  authors  have  used  simplified  versions  to  describe 
liquids  at  high  pressures.  In  the  present  application  the 
theory  is  used  in  its  simplest  form: 


P(V£ 

,T)  =  PR(VS)  +  (W  Eths(Vs)  +  nhkT/V£ 

(4.21) 

EX 

"  Es+Eh 

(4.22) 

Eh 

-  (3/2)nhkT 

(4.23) 

"h 

=  N(V£-VS)/MVS 

(4.24) 

where  subscripts  "s,"  "h,"  "i,"  denote  solid,  holes,  and  liquid 
respectively;  N  is  Avogadro's  number,  M  is  molecular  weight 
and  k  is  Boltzmann's  constant.  Substituting  Eq .  (4.24)  into 
(4.21)  yields 

P(V£,T)  =  Pk(Vs)  +  <rs/Vs>  Echs(Vs>  +  ?,T(VrVs)/MV£Vs  (4.25) 

Equation  (4.25)  does  not  contain  an  interaction  term  between 
holes  and  molecules,  though  that  is  important  in  Eyring's  theory. 
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However  it  may  be  less  important  at  high  pressure  than  proper 
treatment  of  the  solid  compression. 

The  above  equations  are  supplemented  by  an  expression 
for  ,  the  change  of  volume  on  melting  at  constant 

pressure: 

AVm  =  AH(dT/dP)/Tm  =  VrVs  (4.26) 

With  these  equations  the  Hugoniot  can  be  continued  into  the 
liquid  region. 

V.  Melting  Equations 

Attempts  to  predict  melting  parameters  from  atomic 
theories  have  been  many  and  varied.  The  earliest  one  normally 
noted  is  that  of  Lindemann  in  1910  (Ref.  4).  Assuming  an 
Einstein  model  of  a  solid  with  single  vibration  frequency  f, 
suppose  that  the  amplitude  of  vibration  increases  with  temperature, 
and  that  when  the  amplitude  reaches  a  critical  fraction  of  the 
interatomic  distance,  melting  occurs.  Equating  the  energy  of 
vibration  to  the  thermal  energy  of  the  crystal  leads  to  an 
equation  of  the  form 

3R  Tm  =  C  f2  V2^3  M  (5.1) 

where  T^  =  melting  temperature,  =  molar  volume  at  Tm, 

M  is  molecular  weight,  and  C  is  a  constant.  By  writing 

9  2 

4n**f  =  k/m  ,  where  m  is  atomic  mass,  and  setting  the  energy 

2 

of  an  oscillator  equal  to  k(r  -  r  )  /2  ,  where  r  -  rQ  is  the 


deviation  of  interatomic  distance  from  equilibrium,  we  obtain 


e  =  k(r  -  rQ)2/2 
k  =  d2e/dr2  . 


Identify  the  oscillator  energy,  e,  with  the  molar  energy, 
£k,  of  cold  compression  by  the  relation 

£k  »  3Ne  =  £k(V)  . 

Then  Eq.  (5.1)  can  be  converted  to  the  following  relation: 


3M<*2V2(d2Ek/dV2)/2 

(5.2) 

A  V2  d2  Ek/dV2 

(5.3) 

where  V  =  specific  volume  as  in  Section  4,  =  specific 

internal  energy,  »  =  fraction  of  interatomic  distance  at  which 
melting  occurs,  and  M  is  molecular  weight.  Eq.  (5.3)  is  the 
form  given  by  Urlin  and  Ivanov  in  Reference  2.  They  also 
propose  an  alternative  melting  law  in  the  form 


L  V  /T  AV  =  Ra 
mm  m 


(5.4) 


where  L  =  latent  heat  of  fusion,  Vm  =  specific  volume  at 
which  melting  occurs,  T^  =  melting  temperature,  AV^  is  volume 
change  on  melting,  and  Ra  is  a  constant.  This  is  a  modifica¬ 
tion  of  the  Lindemann  formula  which  can  be  seen  as  follows: 

d2Ek/dV2  =  (-dPk/dV)p_0 
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where  the  last  expression  is  obtained  from  Eq.  (4.17). 

Substituting  this  into  Eq.  (5.3)  yields 

RT  =  AV2b,  /V 
I  o 

Differentiate  with  respect  to  P  to  obtain 
RdT/dP  =  2AV(dV/dP)b1/VQ 
or 

dP/dT  =  L/TAV^ 

-  RVo(dP/dV)/2Ab]V  ,  (5.5) 

which  is  of  the  form  Eq.  (5.4)  with  dP/dV  assumed  constant. 

One  of  the  most  commonly  used  forms  for  the  melting 
curve  is  Simon's  equation  (Ref.  5): 

P  -  Pm  +  a  «  a(T/Tm)C  (5.6) 

where  P^  and  represent  one  point  on  the  melting  curve  and 
C  and  a  are  constants  determined  from  the  relations 

aC  =  T  (dP/dT)-  (5.7) 

iu  im 

and 

a  =  <dEk/dV)T=0,P=0  (5.8) 

according  to  Simon.  Later  work  by  Salter,  using  the  Mie- 
Gruneisen  equation  of  state,  identifies  C  as 

C  =  (6rs+i)/(6rg-2) 

where  Ps  is  the  Gruneisen  parameter  of  Eq.  (4.4). 


(5.9) 
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A  melting  relation  recently  proposed  by  G.  C.  Kennedy 
(Ref  6)  relates  volume  on  the  melting  curve  to  temperature: 

Tm  -  Tm°^  +  Ck  W/V  <5.10) 

where  Tm°  is  melting  temperature  at  volume  VQ  and  AV  = 

(V  -Vm) .  Ross  and  Alder  (Ref.  7)  have  criticized  this  as  being 
of  lower  validity  than  Lindemann's  law  and  of  giving  too  low 
values  of  Tm  at  high  compressions.  Gilvarry  (Ref.  8)  suggests 
that  it  is  the  first  term  in  the  expansion  of  the  formula 

T  /T°  =  (V0/Vj2(r o'l/3)  (5.11) 

m  m  o  m 

where  is  the  Gruneisen  parameter  at  V  . 

At  present  it  appears  very  much  as  if  there  is  as  much 
justification  for  picking  one  rule  as  another,  in  the  absence 
of  experimental  data.  In  the  calculations  to  be  reported  later 
we  use  the  Simon  equation  and  the  Kennedy  equation. 

VI.  Calculation  of  the  Hugoniot 

The  computing  process  is  illustrated  in  Fig.  7.  Volume 
and  temperature  are  assumed  to  be  known  on  the  phase  boundary 
and  on  the  Hugoniot  at  pressure  P.  P  is  increased  to  P  +  AP 
and  new  values  of  temperature  are  calculated  from  Eqs.  (5.6) 
and  (4.5).  This  allows  calculation  of  the  coefficients  dV/dP 
from  Eqs.  (1.5)  and  (3.1).  Subscripts  1  and  2  refer  to  values 
on  the  solid-mixed  phase  boundary  and  the  Hugoniot,  respectively. 
Values  of  V^  and  V2  are  then  determined  from  the  relations 


V 


V02 


V0I 


Fig.  7 

Computing  Procedure  for  the  Construction 
of  the  Equilibrium  Hugoniot 
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VX(P  +  AP)  =  VX(P)  +  (.5)((dV1/dP)p+Ap  +  (dV1/ dP)p)AP  (6.1) 

V2(P  +  AP)  =  V2(P)  +  .5((dV2/dP)p+Ap  +  (dV2/dP)pjAP  (6.2) 


The  equations  are  iterated  until  V^(P  +  AP)  and  V2(P  +  AP) 
do  not  change,  then  the  process  is  repeated. 

After  each  set  of  values  (Vp  V2)  has  been  calculated, 
a  test  is  made  to  determine  whether  or  not  an  intersection 
between  V-^(P)  and  V2(P)  has  occured.  If  it  has,  the 
Hugoniot  Eq.  (3.1)  is  replaced  by  Eq.  (3.20)  for  the  mixed 
phase  and  the  process  continues,  testing  at  each  step  to  see  if 
the  Hugoniot  has  entered  the  liquid  phase.  When  it  does,  Eq.  (3.21) 
is  used  again  with  the  equations  of  Section  IVB  for  the  equation 
of  state. 

A  flow  chart  for  the  computing  program  is  shown  in 
Fig.  8;  definition  of  symbols  and  a  program  listing  are  given 
in  Appendix  A. 

The  output  of  the  program  is  illustrated  in  Figs.  9-13 
and  in  Table  3  for  lead.  Fig.  9  shows  the  total  Hugoniot  to 
one  megabar  pressure  using  the  Simon  Equation.  It  enters  the 
mixed  phase  region  from  the  solid  at  about  .645  megabars  and 
leaves  at  about  .675.  Inspection  of  the  curve  shows  that  a 
single  shock  will  be  stable  at  all  pressures  (elastic  waves  are 
ignored).  This  is  verified  by  the  slopes  given  in  Table  3. 

The  narrow  mixed  phase  region  and  the  small  kink  it  produces 
in  the  Hugoniot  curve  suggest  also  that  melting  will  have  a 
small  effect  on  wave  propagation,  even  if  it  occurs  in  shock. 


This  remains  to  be  verified  by  incorporation  of  the  model  into 
a  one-D  wave  program. 


Table  III 


Slopes  of  Hugoniots  and  Rayleigh  Lines 
at  Melting  Phase  Boundaries 


p, 

Megabars 

Solid 

ld?/dVl Hugoniot 
Megabar  g/cc 

Mixed  Phase  Liquid 

|(p-P0)/(v0-v)| 

Megabar  cc/g 

Solid  Liquid 

Boundary  Boundary 

* 

.645 

67.4 

107.2 

21.9 

•  •  • 

.675* 

•  *  • 

114.1  52.6 

•  •  • 

22.7 

.392+ 

39.0 

37.2 

16.0 

*  »  • 

.428+ 

•  •  • 

39.4  42.4 

*  •  i 

16.8 

% . . . — —  —  ~  -  . . . ■  - - 

Simon  Equation  (5.6)  'Kennedy  Equation  (5.10) 

Figure  10  shows  the  Hugoniot  in  the  region  of  mixed 
phase  in  more  detail,  again  for  the  Simon  Equation. 

Figure  11  shows  temperatures  on  the  melting  curve  and 
on  the  Hugoniot  for  the  Simon  Equation.  In  region  A,  where 
the  Hugoniot  is  passing  through  the  mixed  phase,  the  two 
curves  should  coincide.  That  they  do  not  reflects  imperfections 
in  the  equation  of  state. 

Flaws  in  the  equation  of  state  are  revealed  when 
measured  specific  volume  of  the  solid  at  melting,  V  ,  is 
found  to  disagree  with  that  calculated  from  the  equation 

T 

Vol  ■  Vo2  +  iT0l<5V/aT>P=P  dT  >  (S-D 

To2  0 
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where  T  2  is  room  temperature;  VQ2  is  specific  volume  at 
room  temperature  and  pressure  PQ;  Tq^  is  melting  temperature 
at  Pq;  and  is  specific  volume  at  P  ,  Tq-^  .  For  this 

reason  a  procedure  for  calculating  VQ^  is  incorporated  in 
the  program  (Appendix  A).  Even  with  this  correction,  a  slight 
difference  between  TIN  and  T2N  occurs,  as  shown  in  Fig.  11. 

The  difficulty  probably  arises  from  a  minor  inconsistency  between 
the  equations  for  P  and  P^  . 

Parameters  and  material  constants  used  in  the  calcula¬ 
tions  for  lead  are  listed  in  Appendix  A. 

Figures  12,  13  and  14  display  the  results  when  the  Simon 
Equation  is  replaced  by  the  Kennedy  Equation.  The  temperature 
and  pressure  at  which  the  Hugoniot  enter  the  melting  region  are 
about  half  the  values  obtained  with  the  Simon  Equation.  These 
tremendous  differences  represent  the  state  of  our  ignorance  about 
the  melting  process  at  high  pressures,  and  unless  detection  of 
melting  in  shock  is  possible,  that  state  of  affairs  is  very 
likely  to  persist. 


Fig.  8. --Flow  Chart  for  Calculation  of  Hugoniot 
when  Melting  Occurs  ' 


Remarks 


See  Table  A1  for 
symbols  and  defini¬ 
tions  . 


P:  *  P  +  DELP 


Subroutine  PB(C0EF1N 
P,  V1N,  TIN) 


Sub.  RH2  (C0EF2N,  P, 
V2N,  T2N) 


Vln  is  specific  vol¬ 
ume  on  the  solid-mixed 
phase  boundary.  V2n 
is  specific  volume  on 
the  Hugoniot  curve  in 
the  solid.  V3n  is 
specific  volume  on  the 
Hugoniot  curve  in  the 
mixed  phase  region. 

V4n  is  specific  volume 
on  the  mixed  phase- - 
liquid  boundary.  V5n 
is  specific  volume  on 
the  Hugoniot  liquid 
phase.  T  is  tempera¬ 
ture,  °K;  P  is  pressure 
in  megabars. 
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Fig.  S .--Continued 
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Transfer 


Operation 


Transfer 

to 


Remarks 


Control  index  for 
calculating  slopes 


P:  =  P  +  DELP 


Appropriate  subroutines 
(see  listing) . 


transfer 


Operation 


Compute  Sr,  SH  at 
intersection 


sr!  >  ISn^es. 


“Continued 

Transfer 

Remarks 

to 

Sr  =  SLOPER  = 

(p-po)/(v( 
Sb  -  |ap/av5! 

K  =  2 


V5N  ^  V4N  ?  >S2_ 


p  <  pstop  ?  >yes_ 


LIQUID 


another 
scase  to  be 
Xjrun? 


BOUNDARY  BETWEEN  THE 


HUGONIOT  CURVE  FOR  LEAD  IN  THE  MIXED  PHASE  REGION.  SIMON  MELTING  EQUATION 


.50 


IO  o 


to 

fO 


8 


i a 
N 


sy vsvo3W  d  3HDSS3Hd 
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APPENDIX  A 


PROGRAM  FOR  COMPUTING  THE  HUGONIOT  OF 
A  MELTING  SOLID 

Program  Name :  MELT 
Language :  F0RTRAN  IV 

Constants 


Equation 

Program  symbol 

Remarks 

p 

0 

3.1,  5.6 

P0 

atmospheric  pressure, 
megabars 

po 

4.4a 

RH00 

RH02 

density  at  Pq)  T  =  0°K 

density  at  P  ,  T  =  T02 
(room  temperature) 

3R 

4.5 

R3 

R  =  gas  constant, 

Mb  cc/g° 

specific  volume  of  solid  at 

Voi 

A1 

V01 

melting  temperature 
and  P  *  P 

o 

Tol 

'  A1 

T01 

melting  temperature  at 

P  =  P 

o 

aQ >  3^  >  ^2 >  3^ 

4.4 

A0  ,A1  ,A2  ,A3 

bl,b2,b3 

4.3 

B1,B2,B3 

e 

o 

4.9 

TD0 

V 

0 

4.4a 

V0 

AP 

7.1 

DELP 

a  ,c 

5.6 

A ,  C 

Vo2 

3.1 

V02 

sp.  vol.  on  Hugon.  at 

P  =  P 

0 
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Report  symbol  Equation 


Program  symbol 
T02 
PST0P 

IC0NT 


Remarks 


temperature  at  P0,V02 

pressure  at  which  com¬ 
putation  stops 

integer  allowing  several 
data  sets  to  be  used  in 
one  run 


dT  A1  DELT  increment  in  temperature 

used  to  calculate  V01 
from  V02 ,T01 ,T02 


Parameters  In  MAIN 

A  final  letter  "N"  on  a  parameter  symbol  indicates  a  value 
at  pressure  P  +  AP,  an  "0"  indicates  a  value  at  P;  e.g.,  V1N  = 

VI  at  P  +  AP,  V10  -  VI  at  P. 


Report  svrabol 

Equation 

Program  svmbol 

ALG,ALIN 

ASQ,ACUBE 

Remarks 

Program  computed  con¬ 
stants  used  in  evalua¬ 
ting  TD 

V1 

1.5 

V1N,V10 

specific  volume  on  solid- 
mixed  phase  boundary 

T 

1.6,  5.6 

T1N,T10 

melting  temperature 

V 

3.1 

V2N,V20 

specific  volume  on 
solid  phase  Hugoniot 

T 

4.5-4.20 

T2N,T20 

temperature  on  solid 
phase  Hugoniot 

V 

3.18,3.19  . 

V3N,V30 

specific  volume  on  mixed 
phase  Hugoniot 

T 

3.20 

T3N ,T30 

temperature  on  mixed 
phase  Hugoniot 

V2 

1.2 

V4N,V40 

specific  volume  on  liquid- 
mixed  phase  boundary 
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Report  symbol 

Equation 

Program  symbol 

Remarks 

T4N,T40 

same  as  T1N,T10 

V. 

z 

4.25 

V5N,V50 

specific  volume  on 
liquid  phase  Hugoniot 

T 

4.25 

T5N ,T50 

temperature  on  liquid 
phase  Hugoniot 

dVj/dP 

1.5 

C0EF1 

slope  of  phase  boundary 

dV/dP 

3.2 

C0EF2 

slope  of  solid  phase 
Hugoniot 

V 

A1 

VX 

value  of  V  at  P  =  P  , 

T  =  TX  (T02^TX^T01)° 
used  to  calculate  V01 
from  V02 ,T01 ,T02 

T 

A1 

TX 

temperature.  T02^TX^T01 
used  in  the  calculation 
of  V01 . 

TXT 

TXT  =  TX-DELT 

VXT 

value  of  VX  at  TXT 

MU,TD1,X,D,DP  , 
GI,Q1,W,ET,M1,  ) 
ETT ,R1 ,S1 

same  as  in  subroutine  PB 

JC 

integer  used  to  control 
what  part  of  a  subroutine 
is  to  be  used  in  a  partic¬ 
ular  calculation 

K 

controlling  integer  used 
in  main,  similar  to  JC 
in  operation 

-dP/dV 

3.1 

SLOPEB 

dP/dV  for  the  Hugoniot 
just  inside  the  one  phase 

~dP/dV  3 

.20,3.21 

SLOPEH 

dP/dV  for  the  Hugoniot 
just  inside  the  next  phase 

3.3 

SLOPER 

slope  of  the  Rayleigh  line 
from  the  foot  of  the 
Hugoniot  to  its  intersec¬ 
tion  with  the  phase  boun¬ 
dary 

1  , 

E:  &  I 


l 

ft  -t  Of 

f 

}l  * 

S  ‘l 
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s> 

{ 

•*/' 

l 

Report  svmbol 

Equation 

Program  svmbol 

Remarks 

1 

y. , 

dV/dP 

3.20 

C0EF3 

slope  of  mixed  phase 

1 

r 

tv<<' 

**? 

A 

Hugoniot 

£ 

V  ' 

dV/'dP 

3.21 

C0EF5 

slope  of  liquid  phase 

’ 

Hugoniot 

i 

1 

f 

yO 

VI 

sp.  vol,  pressure  and 

5 

temperature  at  inter¬ 

l 

PI 

* 

section  of  V-.  (P)  and 

! 

jt 

A 

S/,(p)  1 

•vVn 

TI 

<u 

1 

=> 

S? 

VI2 

sp.  vol.,  pressure  and 

j£ 

j 

, temperature  at  inter¬ 

? 

' 

PI2 

* 

section  of  V0(P)  and 

'l 

f* 

V(p)  3 

K 

5 

* 

TI2 

H* 

M 

_ / 

| 

dP/dT 

1.1 

PT 

Clausius-Clapeyron 

i 

coefficient 

l 

AV 

1.1,  1.2 

DELV 

sp.  vol.  change  on 

X 

f 

melting 

£ 

fe 

AH 

1.1 

LATHT 

latent  heat  of  melting 

Parameters  in  SUBROUTINE  PB(CCfEFlN.P.VlH.TlN) 


Report  svmbol 

Equation 

Program  svmbol 

Remarks  j 

V1NN 

jj 

temporary  value  of  V1N  '  ' 
used  to  test  convergence  | 
of  iteration  ? 

4.4a 

MU 

l 

e 

4.9 

TD1 

i 

1 

T 

5.6 

Tl 

( 

X 

6/T 

dT/dP 

CC 

4 

reciprocal  of  Clausius-  | 

Clapeyron  coeff J.cient  | 

D 

4.7 

D 

1 

Debye  function  1 

D’ 

4.16 

D? 

j 
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Report  symbol 

Eauation 

Program  symbol 

Remarks 

r 

4.4 

G1 

cv 

4.12 

CV1 

d(r/v)/dv 

4.13 

Q1 

dPk/dV 

4.17 

W 

Eth 

4.5 

ET 

(*p/°v)t 

4.18 

Ml 

Cv  4 

.10-4.12 

ETT 

<*F‘th/dV>T 

4.15 

ETV1 

(3V/SP)t 

•  •  • 

R1 

(9V/ST)p 

4.19 

SI 

dVx/dP 

1.5 

C01 

same  as  C0EF1N 

G010 

same  as  C0EF10 

Parameters  in 

SUBROUTINE  RH2(V.P.CDEF2N,T2N) 

Report  symbol 

Equation 

Program  symbol 

Remarks 

V 

V2NN 

temporary  value  of  V2N 
used  to  test  convergence 
of  iteration 

u. 

4.49 

MU 

V 

3.1 

V 

same  as  V2N 

r 

4.4 

G2 

e 

4.9 

TD2 

pk 

4.3 

PK0 

D 

4.7 

D 

Gv 

4.12 

CV 

Eth 

4.5 

ETH 

dPk/dV 

4.17 

PPK 

b  -f  r- 


B  ' 


?r  f 


t  ! 
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Report  svrabol 

Equation 

Program  svmbol 

Remarks 

d(r/v)/dv 

4.13 

Q2 

(3Eth/SV)T 

4.15 

ETV2 

(3P/av)T 

4.18 

M2 

(sp/sv)s 

4.20 

MS2 

dV/dP 

3.2 

C0EF2N 

e/T 

4.5 

X 

D' 

4.16 

DX 

Parameters  in  SUBROUTINE  RH3(rriRim  p  wsm  v4N^ 

Report  svmbol 

Equation 

Program  svmbol 

Remarks 

V3NN 

temporary  value  of  V3N 
used  to  test  convergence 
of  iteration 

T1 

temperature  at  P  +  AP, 
transferred  from  SUB¬ 
ROUTINE  PB 

VI 

V1N  at  P  +  Ap  transferred 
from  Sub.  PB 

d2P/dT2 

3.19 

PTT 

A 

3.19 

Y1 

•5(P-P0)- 

TdP/dT 

3.20 

ANUM1 

A-B 

3.20 

DENI 

C0EF3N 

same  as  in  MAIN 

V4N 

n  ii  it  u 

C0EF30 

ti  ii  it  ii 

.1 


i 


1 
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Parameters  in  SUBROUTINE  TEMPl(P) 

Report  symbol 

Equation  Program  symbol 

Remarks 

dP/dT 

5.6  PT 

d2P/dT2 

5.6  PTT 

T 

5.6  T 

T  =  TIN  of  the  main 
program 

Values  of  constants  used  for  lead: 

po 

=  11.616  ;  3R 

0.1202506  x  10‘5 

1! 

r-» 

o 

H 

601  °K 

Vo2 

=  0.08818340  ;  Tq2  = 

293°K 

a 

0 

=  2.7091  ;  a1 

-2.5282 

a2 

=  1.413  ;  a^  = 

0.0 

bi 

=  0.54168  ;  b2  -  0.749041 

;  b3  =  0.605839 

i> 

O 

=  96.3 

a 

=  0.06257  c 

1.20436 

Constants 

independent  of  material: 

Po 

-  6 

=  1.034  x  10  megabars 

AP 

=  .005  megabars 

PST0P 

=  1  megabar 

Computation  of  VQ^ 

The  equation  of  state  used  for  the  solid  in  Section  IVA 
is  self-consistent  but  not  entirely  consistent  with  all  available 
measured  data.  In  particular  if  V0^  is  taken  to  be  the  handbook 
value,  t  T2  at  the  intersection  of  V-^P)  and  V2(P) .  To 
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remedy  this,  VQ^  was  calculated  from  the  equation 

T 

Vol  -  Vo2  +  jT  ('6V/aT)P*P  dT  (A1) 

To2  ° 

where  T02  is  the  value  at  the  foot  of  the  Hugoniot  and  TQ^  is 
the  handbook  value  for  melting  at  P  *  PQ  .  Handbook  and  calcu¬ 
lated  values  of  VQ^  are: 

VQ^  ®  .091148  cc/g  (metal  reference  hbk) 

Vol  *  .090352  cc/g  (Eq.  (Al)) 


Program  listing 

The  program  listing  follows  on  separate  pages. 
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//LISTJCBS  JC2  (0000, CC, 5,  10  J ,  PACKRCC.#  ,PSGLEVEL*l 

//JCBLI  B  0C  UNIT»231  l*VGLUMfc»$ER*DLIB02,0IS?aOLD,DSNA‘'E*$YSl.wTlL.TY 
//STEP  EXEC  5GP=L  ISTCRD 

//SVSLST  00  SYSC*UT=A  *  CCfi  -  ( LRGCLa80, BLKSl 7  S  '80»RECFM=F  1 
//SYS004  DC  L'M  IT=SYSOA,  VCLUPE=SER*SCR001..C*NAME=L  t  STGR, 

//  SPACE=(CYL,(5,1) 1 ,0C0= ( RECFM=  FBS, OLKSI ic  =30 , LRECL=80) 

//SY.SRCR  CO  * 


(6F236I 

ALLOC. 

FOR 

LISTJOBS  STEP 

IEF237I 

JOHLIB 

OS 

293 

IFF 237 1 

SYS004 

ON 

290 

IEF237I 

SYSPCR 

0?I 

OOC 

OUVALL 


PAGE  1 


THIS  PROGRAM  INTEGRATES  THE  EQUILIBRIUM  HUGONICT  P-V  CURVE 
CF  A  SOLID  IN  INITIAL  STATE  ffi2iP02(*PC) ,  VQ2  ASSUMING  A  SINGLE  SHOCK 
FROM  PO  TO  THE  FINAL  PRESSURE  P.  IT  SIMULTANEOUSLY  CALCULATES 
THE  MELTING  CURVE  AMO  AT  EACH  STEP  TESTS  TO  SEE  WHETHER  MELTING  OCCURS. 

IF  IT  COES»  THE  SLOPE  OF  THE  HUGONIOT  IN  THE  MIXED  PHASE  IS  COMPARED 
WITH  THE  SLOPE  OP  THE  RAYLEIGH  LINE  TO  SEE  WHETHER  OR  NCT  A  SINGLE 
SHOCK  IS  STILL  STABLE.  IF  NOT,  THE  COMPUTATION  IS  TERMINATED: 

IF  IT  IS,  THE  COMPUTATION  IS  CONTINUEO  THROUGH  THE  MIXED  PHASE  REGION 

ANC  TESTS  FOR  INTERSECTION  WITH  THE  BOUNDARY  BETWEEN  THE  -«IXED  PHASE 

ANC  THE  LIQUIO  PHASE  ARE  MADE.  IF  THE  INTERSECTION  OCCURS  THE  STABILITY  TEST 

IS  REPEATED  AMO  THE  COMPUTATION  IS  STOPPED  IF  INSTABILITY  IS  INDICATED! 

OTHERWISE  IT  CONTINUES  IN  THE  LIQUID  PHASE.  IF  THE  CURVE  RE-ENTERS 

THE  MIXED  PHASE,  THE  COMPUTATION  IS  STOPPED;  OTHERWISE  IT  CONTINUES 

UNTIL  P*PSTCP. 


TO  PUN  THE  PROGRAM,  PREPARE  OATA  ACCORDING  TO  THE  FORMAT  IN 
STATEMENT  NUMBE°S  130  ANC  101 (MAIM) , 

PO- INITIAL  PRESSURE*  1.0346-06 
RHGO*OENS l TY  AT  e0,Z-?0  CFGREES  KELVIN 
rhc2=censity  at  the  nor  cf  the  hugoniot 
MOLWT*MCLECULAR  WEIGHT  OF  THE  MATERIAL  IN  GRAMS 
T02= TEMPERATURE  AT  THE  ':COT  OF  THE  HUGONIOT 
TCC=CEBYF  TEMPERATURE  <>  •  VO*1/RHOO 

AQ.Al.A2, A3  ARE  THE  COE.  FICIENTS  IN  THE  EQUATION  FOR  GRlJVEISEN'S  PARAMETER,  G 
WHERF  G=  AC+MU**  C  A l «  MU* (  A 2 ♦MIJ* A3  ) )  AND  MU=VC/V-1. 

ei.B2.B3  ARE  THE  CCEFFIC  IcMTS  IN  THE  F.OUATICN  FOR  PRESSURE  ON  THE 
ZERO  CEGREG  ISOTHERM J 
PK=MU<'(0UMU*(I32+MU-*B3) ) 

OFLP  IS  THE  INCPEMFNT  IN  P  USED  FOR  THE  I NT6G.R \T ION. 

ICPNT  IS  AN  INTEGER  WHICH  ENABLES  THE  USE  T  TU  INPUT  AS  MANY  DATA  SETS 
AS  QESIREC  IN  ONE  RUN  OF  THE  PROGRAM, 

IC0NT*1  OR  2  FOR  A  CC^PLETE  OR  PARTIAL  NEW  SIT  Or  DATA  (RESP.)  TO  BE  READ  IN. 
IN  THE  LAST  DATA  SET  IC0AT*AN  INTEGER  OTHER  THAN  1  OR  2  FOR  THE 
FRCGRAM  TO  TERMINATE  AFTER  THE  LAST  EXECUTION. 

VPl,  THE  SPECIFIC  VOLUME  OF  THC  SOLID  AT  TCI, PC  IS  CALCULATE!) 

IN  THE  PROGRAM  SO  THAT  TEMPERATURE  ON  THE  MELTING  CURVF  IS  COMPATIBLE 
WITH  THAT  ON  THI  MUGCNIOT. 

COMMON  AQ, A1 , A2 » A3 , B1 »D2, B3, ALG, AL IN.ASQ, ACUBF , VO »P0 »DGLP * R3 » VIC, 

+V20,VC<  , :02rTD0 ,T0l,T20,JC,PTT 
CCMMON  /a’RT  1/  P , V IN, T IN ,CCEF  1N.V2N.T2N.CCEF2N 

COMMCN/CGIPSi/TDl , G 1 , CVl , 01 , V, ET ,M1 ,ETT ,ETVI , Rl , SI , PT , MU1 ,Xi,D!, 
hOP,TT,VV 

CCMMCN/Cn2?iF2/r.2f  TI)2,CV,ETH,Df'K,Q2,FTV2,M2,MS2,MIJ2tX2,02 
COMMON /PRESOK/PKO 

REAL  L  4THT ,  MUl ,  MIJ2  ,  M 1  ,M2,  MS2 ,  MU,  MCLWT 
WRITE  16,10) 

10  FORMAT  PI1/'  •,8X,'P,,13X,)Vl»,lSX,'Tl',1.3X,'V2',13X,'T2»i 

103  R  6  AC  1 5, 100)  PO,  pMCO  ,RH()2,  MOLW'T ,  T01 ,  T02  »  TDD 

100  FCPMAT(E20.7,fcF10.6) 

104  R  E AO  I S , 101 )  AO, A l , A2, A3, 9 1 , B2 , B3,0ELP, PSTOP, I CONT 

101  FORMAT! 7F10.6/2F1C. 6, 12) 

V0=1./RHC0 

VC?*l./RH02 


51 


DUVAL'. 


1U=('S,955*'.1°F-51^,GIMT 

iLG»AC-Al*A2-A3 

ALJN=Al-«2«-A3 

ASC*(A2-A?)/2. 

ACUBE=A3/3. 

0ELT*0.5 

VX*VQ2 

TX*T02 

<56  PU=VC/VX-1. 

TCl*TCr-1U) 

X-TD1/TX 

D=l.-X*(0.37B-0.C5*Xl 
CP=-0. 375+0. l*X 
C-l  =  G(Pl>'» 

C1=G<VX,V0)  ^ 

w— _ ( VG - (  2  > )  /( VX**2> 

FT=R3*C*TX 

t.'l*s+ST*01-R»«iTiH*JI  (Gi/vxJ^moP 
0TT*R3*<  1.-0.C5*!  >**2M 
Rl=l./’U 

Sl=-fiHGl*ETT/VX 

VXT=VX 

TXT=TX 

VX=VX*-$1*DEL  T 
TX=TX  tCCLT 

UMTX.LG.TCH  GC  TO  9€ 

7  ppUIv  { e20.rt'»l»0S620,7.'-RH0C»fE20.7,«»«3», 620.7,  »-V01«// 
4T20.7, 'sTCl' .G20.7, •aAO*t62C.7»,«Al,»G2rt*7**»A2'// 

+  C20.  7,  «=  \3' ,l:?0  .  / 1 '  =>D  1 ' » F20,  7i 'aB?' iE?0.7,  =  IH  // 

*r.20Ht  ,®TD0,»f?0.7,  ‘»CSL?'*e?.0.7t  ,ST02>) 

+  ,l:/'’0. 7»  '=VO'  »E2C. 7» '  =\02'  ) 

v;n--vc2 

tu;=tci 

T  2U  =  T02 
P  =  I>C 

CALL  TEV?  1 1  P ) 

CC£F2N=“V01/B1 

JO  1 

V1A-VCI 

CALL  TG  !P ( P i T IN » PT  ) 

CALL  ?•’  (  V  IN  i  PC  *  C-3  rf  IN  »T]N  ) 

L’.*L'  rh?.(V2N«Pr  ,CrEp2N»T2Nl 

\RHHU:,244  7)>,tIN»FT  ,VlN  iMl'l  *  T01 .  XI  ,<U  »0P»r.l  »CV1  *01 .  FT  i  ’-'l »  ETT » 

l‘,T°  t  TG  2p4 »  T2*!»  V <ri;ML'2  »T02  » X2  «0 2 1  CPK»G2 » CV*  03 » E TH  i  v<?  *  F.T2V2  • 

i'  MS2,Ci:EG2N,P<n 
t,  VIC  =  V1.\ 

V 20= VPN 
T  lr  =  T  IN 
T  2C=T2\‘ 


i 
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5  P-P+MLP 
CALL  TEMP(P,TIN,PT) 

CALL  PG(V1N»P»CG EFIM, TIN) 

CALL  PP2(V2N,P,C0EF2NtT2N) 

WRITEI6,2447)P,TIN,PT,V1N,MU),TD1  »Xl,D5  ,DP,G1  ,CVl  ,*31  ,ET,H1 ,ETT, 
l  ETV1,R1,SI,  CGEF1N 

WRITE (6, 2448)  <> ,  T2N,  V2N.MU2 , TD2,  X2 , 02,  DPK, G2 ,  CV,  Q2 , ETH.M2  ,ET2V2 , 
l  PS2,CCEF2N,PKO 
C 

C  TEST-  FOR  I NTERSECT.I ON  WITH  PHASE  3CUADARY 

C 

IF  IV2N  .IT.  ViN  )  GC  TM 
C 

C  FINO  INTERSECTION  OF  R-H  AND  PHASE  3CUNQARY 

C 

V  I- {  V2N*V10-VlN*'/20)/(V2N-V20-VlN+Vl0J 
PI=»P«-CELP*(  VI-V1N  )/(  V1N-V10) 

Tl*TlNMTlN-Tin)*m-VlNj/(  V1N-V10) 

WRITE  (6,16)  PI, TI, VI 

16  FORMAT  I'0',G20.7,'*Pt‘,E20.7,'’TI',E20.7,'=Vl'J 
SLOPE 6SCELP/ ( V20-V2N ) 

V»R  l T E ( 6 , 102  ) 

102  F  CP  HAT  ( '  C*  /»*/'  *,20X,  'CONTINUE  IN  MIXED  PHASE'///) 

C 

C  CONTINUE  IN  MIXOO  PHASE 

c 

p=p  I 
VIN»VI 
V3N=V I 
T1N-TI 
JC-1 

CALL  TEMP(P,TIN,PT) 

CALL  PfMV'N.P, COFFIN,  TIN* 

CALL  PF3  (V3.N,  P  ,  COFFIN  ,  VAN  ) 

WRITE  (6, 2447 IP,  T !  N,  PT , VIN ,MUi ,T01 , XI ,0 1 ,0P ,G1  ,CVl  ,-11 , ET ,M1 , ETT , 

1  ETV1  ,R1  ,S1  ,  COEMM 

WRITE  (6,17)  P, T  IN, V  1N,C0GF1N,C0EF3N, V4N 

17  FORMAT  (F20.7, ',S2C.7, 'aTlN* ,E2C.7, • » VIM • /E20. 7, ' *CCEF IN ' , F20. 7 
♦  ,  •»CG6F3N',6?0.7,  '*•/ 4N'  > 

K=1 
JC*2 

11  V10=V1N 

■  ■  .uni 

V  3U*»  £i\ 

V40*V4i\ 

T 1C=T IN 
P*P+CELP 

CALL  T  EMP (P , TIN, »T ) 

CALL  P0 I V1N , P, COEF IN, T IN ) 

CALL  RH3  I V3N ,  P ,  C'JEF  3-N  ,V4M) 

WRITE(5,2447)P,TiK,OT,ViN,MUi,T0l,Xl,01,DP,Gl,Cvi ,31 ,ET , HI , ETT , 

1  ETV I , R 1 , S 1 ,  COFFIN 

W°!TF  (6,13)  0 , V  IN , V  2N , V4N, TIN 

18  FORMAT  IE20.7,  •=?•  ,F20.7,  '=»V1N • , E20.7,  * -V3N* /E20. 7 , ' »V4N»  ,  E20.7 ,  •  a 
♦TIN' ) 

IFfKvKE.il  GO  TO  55 
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SLCPEH=CQLP/(Vt-V3N» 

SLf.'PEP=lFl-P0}/(VC2-V!{ 

C 

C  SLOPES  IS  THE  SLO?C  f.F  RAYLEIGH  LINE  FP.CM  THE  PCCT  OF  THE  HUGCMOT 
C  TO  US  INTERSECTION  WITH  THE  PHASE  EDUNOA'RY 

C  SLCPEH  IS  OP/DV  FOP  THE  HUOONIOT  JUST  INSIDE  THE  MIXED  PHASE  REGION 
C  SLCPEH  IS  DP/OV  FOR  THE  HIGONIOT  IN  THE  SOLID 
C  JUST  OUTSIDE  THE  MIXED  PHASE  REGION 
C 

WRITE  I  6, 51)  SLOPCH, SLCPER  ,  SLOPEB 

31  FORMAT (' O* ,615. 7, '=SLCPE.I',5X,E15.?,'=SLCPER' ♦5X.G15. 7, '’SLOPEB •/) 
IF( SLCPEH. GE. SLOPES)  GG  TO  55 
RETURN 
55  K=  2 

IF  l V3N  .LT.  VAN )  GC  TO  U 

V 1 2  =  (VA0*V3N  -  V?0*VA'J)/C  V3N  -  V3C  -  VAN  ♦  VAO) 

PI?  =  P  +  OGLP* (VI?  -  VAN ) / ( VAN  -  VAO) 

T  1 2  *  TIN  +  (TIN  -  T 1C)*( VI2  -  V4N)/(V4N  -  VAC) 
l/PITE  (6,19)  PI2/VI2,  rtZ 

19  FORMAT  (F.20. 7,  '=■>!?•,  E20.  ?,  '-V.?*  ,F?0.  7,  •  -  TI2  ' ) 
l.'  R  I T  p  ( 6  *  1C6  ) 

106  FORMAT  I'C'///'  ',  ?0X, 'CONTINUE.  IN  LIOUJO  PMASIT'///) 

SLOP CG  =  0£LP/ C  V30-V3N') 

V5N  =  V  I  2 
T 5 N’T  12 
P  =  P  I  2 
TAN=T  12 
T IN’ TAN 
VAN= V  I  2 
T  2N=  T  5N 
JC  =  1 

CALL  T  c  M  P ( P , T l N  <  P  T ) 

PTLV’LATHT  (T  IN  ,!>)/<  T1N*PT  * 

V ia’VAN-CELV 

CALL  FEIV1N,P,C<'GFIN,T1N) 

V2N=V5N-CFLV 

CALL  PH? (V?N, P , COEF 2N , T2N ) 

HR  IT  F  (  6 , 602  )OEt.V 
602  FORMA  T ( •  »,filA.6,'=DFLV') 

WRITE(6,2AA7)P, TIN, PT,V1N,MUI,TD1,X1 ,01 , DP, 51 ,C VI ,01 ,ET  »Ml ,ETT , 
l  ETVl  ,R 1 , SI »  COtFIN 

WR I TE (6, 2AAd )  P , T2N, V2M.MU2 ,TD2 »X2,D2»DPK,G2,CV,02,FTH,M2»ET2V2, 

1  MS? , CCEF2N , PKO 
K=1 
JC  =  2 

25  T20-T2N 
T  U  =  T  1 N 
VIC’VIN 
V2C=V2N 
P=P+CI IP 

CALL  TF“P(P, TIN.PT) 

CALL  PEIVIN, P, COFFIN, Tl”> 
l)ELV  =  LATHT(m,P)/lTlN"PT) 

VAN=V 1N+CELV 
TAN’TIN 
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CALL  P-F2iV2N*P,  CQEF2N,T2NI 

V5N«V2N>CELV 

Tii\uT2N 

IF(K.NE.l)  60  TO  70 
SLCPEP«{PI2-P0)/!Vn2-VI2J 
SLOPEMDELP/I  VI  2-V5N) 

WRITEI6,51 I  SL0P6K, SLOPE* »SLOPEB 
I  Ft  SLCPEF.GC .SLOPER )  CO  TO  70 
RETURN 
70  K»2 

V!R!TEt4,244  7)Pf  fiN,PT,VlN,MUi»T01,Xl,0l,0P,Gl,CVl,Gl,ET»Ml,ETT, 

1  ETVl.Rl.Sl,  C0EF1N 

WR I TE ( 6,  2449  >  P  ,  T2N,  V 2K> .  MU2 , T02 .  X2 , 02 ,  DPK.G2  »CV  ,02  ,FTH ,M2 , ET2V2 , 

1  MS2*COEF2N#PKO 

2447  FORMAT  I » O'  ,615.7,  •  *P»  ,5X ,  F15*  1,  •‘•TIN* ,  5X,  Cl 5. 7,  •  =  PT*  , 5X, E15  .7, '  *V 
UN',5X,G15.7,'*MU1'/*  ',E15.7,  '«T01  •  ,5X  ,E  15. 7  , '  =»Xl  • ,  5X,G15. 7,  •  =01  • 
2, 5X, FI  5. 7, '-CP* ,5X,E15.7, ' » G15 , 7, '*C VI  * , 5X , E 1 5 .7 , » =01 • , 5X , 
3E15.7,  '»GT',5X,G15.7,  •■n*.5X,EI5*7»,»ETT»/»  ',E15.7,  •>ETVl'*5X>Et 
45.7,  »=Rl'  ,5X,F15.‘7,  •=  SI ' ,  5X.E1 5. 7, •  <*OOEF1N' ) 

244fi  FGRMAT('C',E15.7,  »  =  P ' , ?X.  E15,  7:  •  »T2NM  » 5.X ,E15.  7,  •  *V2N'  ,SX,E15.7, »=M 
1U?',5X,E15.7,'  =  T:>2'/'  ' ,F IS. 7, *sX2* ,5X,E15. 7,  '*32 ' ,5X ,515.7 ,' =OPK* 
2, 5X,E15.7,'*G2* ,5X,G15.7, *=CV'/'  • , EI5 . 7 , « *02 '% 5X ,Gl 5 . 7, • «GTH« ,5X, 
3E15.7, '-M2»»5X.ei5.7, '»6T2V2' « 5XtE15. 7, • *KS2 • /•  '  ,E1 5.7  , ' =CCEF2N" , 
45X,£15.7,'=PK' J 

WR I T E  t  6, 600 ) P *  V IN » V4N  ,V2N ,V5N, T4N, T5N 
600  FORMAT  (•  '.7EI5.6J 

IFIV5K.GE.V4N)  G3  TC  T> 

RETURN 
72  CONTINUE 

IFtP.LT.PSTUPJGC  TO  2* 

IF(  ICONT.EQ.I)  GO  TO  103 
IFUCCNT.EQ.2)  GO  TO  104 
RETURN 
ENO 

SUBROUTINE  TEMPltP) 

COMMON  AH,  A 1 »  A2,  A3,  B1  ,B2,  B3,  ALG,  ALIN»A S0»  AC'JBf ,  VO  ,P0 , DEI.P?  R3 .  VI 0, 

♦  V2P,  VG2,  TI12 » TOO ,  TCI , T20,  JC»PTT 

C  A  ANC  C  ARE  THE  CONSTANTS  IN  THF  SIMON  EQUATION  OF  MELTING! 

A-0, 06257 
C»l.  20436 
W  R I T  E ( 6 , l)A,C 

1  FCPMATt*  * , 615.7, *aA*  , 5X,015.7* '*C' ) 

RETURN 

ENTRY  TEMPt  °,T , PT  ) 

IF(JC.NF.l)  GO  TO  2 
PT»A*CMT**(C-1.)  )/(TCl**C) 
l>TT»C*(C-l.  )#(P-»C«-A)m**2) 

RETURN 

2  CONTINUE 

T  =  TCl*t(  (P-°0*A)/A)**(l,/Cn 
PT=A*C*tT**<C-i.)  )/(TCi**C) 

PTT=C*(C-1 .  )<MP-0C  +  A  )/{T**2) 

RETURN 

END 

SUBROUTINE  PB( VI, P.C01 ,Tl  1 
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COHCN  AC, Al, A2,  43, ei,B2,B3,ALG«AL!N,ASO, ACU3E,VC,Pn,DELP,R3,VlC, 
♦V20, VC2 , TG2  »  TDC»TC1 » T20, JC,PTT 

CC^MCN/CClP31/T0l,Gl,CVltCl,W,ET,Kl,ETT,ETVl,Pl,Sl,Pr,^Ul,Xl,Dl, 

+CP,TT,VV 

REAL  PL',-U,HU4,MU1 

ifcjc.co.h  coi=o. 

CC10=CC1 
CC  8  J=l,10 
IFCJC.E0.il  GO  TO  10 
V  INN3 V l 

V1=V1C+0.5*(C010+C01 )*D£LP 

1 F( A8SC  C VI— VIMN1/V1 1 .GE.A.E-51  GO  TO  10 

RETURN 

10  .ML=(VC/V1)-1  . 

TC1=TC(MU) 

X=TC1/T1 

CC=1./PT 

0*1. -X*(0. 375-0. OSOXJ 
0P=-0. 275+0. 1*X 
G  1  =  G  ( -VU  I 

CV1  =  R3*(/,.<!C-3.»X/(EXPCX)-1.  n 
01=0 (VI, VC) 

W  =  R  1  +yu<‘  C  2.*32+3.*83*PU) 

W=-VG*W(V1**2  ) 

ET=P3*OTl 

f'l=W+ET*Ql-R3<-TCl*(  (Gl/Vl  )**2)*0P 
FTT  =  R3M  l.-O.OSM  X**2) ) 

ETVl  =  R3*CP-M  C-G1*T01)/V1) 

R 1=  L . /Ml 

S 1  =  -R 1 *G1*  ETT/V 1 

CC1=RI+S1*CC 

TT  =  Tl 

VV=V1 

KL1=MU 

X1  =  X 

01  =  0 

(F(JC.EQ.l)  GO  TO  12 
8  CONTINUE 

WRITE (6,600)  V1.V1NN 

600  FCPMAT ( 'O' , ' SUOR.P0,  ITERATION  F A 1 LEC ' , 2E20. 7 ) 

CAtL  EXIT 
12  CONTINUE 
RETURN 
ENC 

SL9RCLT INE  RF2(V,P,COEF2N,T2N) 

CCRMCN'  AC,A1,  A?,  43,61  ,3  2 , 33 , ALG, AL IN, A  SO, ACU3F , VI »  PO , DELP , R3  ,V10, 
+V20,Vr2,TC2,TGG, TCI , T2C.JC.PT T 
C07*HCN/CC2RF2/G2,  T02, CV , E TH, CP* ,02 , E TV2  ,M2,’1S2,  VU>  ,X2 ,32 
CCMVCN/PRFSUK/PKO 
REAL  »‘L,.M2,mS2,PU2 
l c ( JC . FO , 1 )  C0EF?N=0. 

CCFE2C=CCCF2N 
DO  8  J=l,10 
iF(JC.EQ.l)  GO  TO  10 
V2NN=V 
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V«V2C4C.5*IC0EF2V-»C0EF20)*fiELP 

I P ( AIS ( ( V-V2NN )  /  V )  .G£  .A»E-5J  GO  TO  10 

RETUPfS 

10  MUVO/V-1. 

G2*G(MU> 

TC2*TC(MU) 

PKC*PK(MU) 

ETM(IP-FKC)*V)/G2 

X=nETF/(F3*rC2)+-.3  75)-SQt!Tt  (F.TH/(R3*T02)  +  .3?5J-x*2-0,2)  )/0.l 
T2N-TG2/X 

D=l,-C.375<>X40.C5*X**2 
0CCX=-0. 27540. l*X 
GTV2=P3«G2*TP2*(4-C.27*-.1*X)  /V 
CV=R3* ( 1 •-O.C5’M  X**2) ) 

DPK=B  1  *MU*  <2.*P24VU*3  .*133) 

DPX=-CPK*Ve/(V**2) 

Q2*Q( V  i  VC ) 

H2»CPK4ETF*C2+G2*rrV2/V 
MS2-'M2-CV*T2M*<  <G2/V)**2) 

CCEF2N=  (1«-G2*(V02-V)/(2.*V))/(MS2+G?*{P~PU)/(2.*V)) 

MU2*MU 
X2*X 
C2  =  C 

IF(JC.EO.l)  GO  TO  12 
8  CONTINUE 

WRITE  (6 1 600 ) V*  V2\'N 

600  FCPMAT  ( *  0  * »  'SIJPP  .RH2*  ITERATION  FA  I  LEO*  ,  2E20. 7 ) 

CALL  EXIT 
12  CONTINUE 
RETURN 
END 

SLORGUT l NE  RIO  ( V3N , P , COSF  3N , VAN  ) 

COMMON  AC,Al,A2,A3,f»l.B2,  S3 ,  ALG,  AL  tN,4 SO, ACUOC, VCt  PO, 05 IP , R 3 , V 10, 
+V20,VC2,T02,Tn?,TCl,T20,JC,PTT 

CCNMCN/CC1PI31/TC1,GI,CVI,01,W,ET,mi,FTT,ETVI ,P1 ,Sl ,PT,«Ul ,Xl,01, 
+DP,TT,VV 

PEAL  WU » LATH  t  u  1 1  BU  l 
IF(JC.EO.l)  CnEF3N'=C. 

CCEF3C=CCFF3N 
V3C=V3N 
CC  8  J=l,10 
IF(JC.EQ.l)  GO  TO  10 
I  V2N’N=*V3N 

V2N-V30  ♦  ,5*<CCEF3C  ♦  C0CF3N)*DELP 
I F { AfiSI 5 V3N-V3NN ) /V3N  ).SE • A. E-5)  GO  TO  10 
RETURN 
10  T  1-=TT 
V1*VV 

Y 1 = ( CVl M  1.  ♦'31*T1*$1/V1)  +  T1«<  V3N-V1  )*°TT)/I>T 
DENl=Yl-2.*Tl*Sl-Tl*Rl*PT-IV02-V3N)/2. 

ANUM1».5*{P-P0)-T1*PT 
CCEF3N*0EN1/AMUM1 
VAN=V1+LATHT(T1,P)/(T1*PTI 
IF(JC.EO.l)  GO  TO  12 
8  CONTINUE 
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WRITE (6,  £00) V3N» V3NN 

600  FORMAT  t  *0! , * N5-3-  ITER  .  FAILED' ,  2E2G,:  7) 

CALL  EXIT 
12  CONTINUE 
RETURN 
END 

FLNCTICN  G(MIJ) 

COMMON  AO.Al, A2,A?,81 , B2 , P3 , ALG , AL IN » A SQ , ACIJBE, VO , PC, DELP , R3 tVIO, 

+  V2P*  VC2 «  T02 »  TOO , TCI »T  20, JC, PTT 
REAL  MU 

G  =  AO*vU*(AH’‘UMA2  +  mu*A3)  ) 

RETURN 

END 

FUNCTION  TO (MU) 

CCVNCN  AQ,A1, A2»A?,81,B2.B3,ALG,ALIN,ASO,AClJ3E,VOfPn,OELP.R3,ViO* 
*-V2f,V02,  T02  ,  TOO,  TCI  ,  T2G,  JC.PTT 
REAL  MU 

TC=TOC*EXP(  ALG<-ALCG‘lML*i.  >  *«!.'*(  AHNfMU*(  A  SQ+ML*  AC U»E  )  )  ) 

RETURN 

END 

FUNCTION  C(V,VX) 

COMMON  AC,  Al, A2, 43,31,32,P3»ALG;Al.IN,AS0, ACU8F,  VO , PO , CELP ,R3 , VIC, 
+V2C.V02, T02, TOO, TCI ,T 20, JC.PTT 
REAL  VU 
ML'=IVC/V)-1. 

0=-(V0/(  V<=*3)  )*<A]+MU*(2.+A2:+3.*A3*MU)  )-G(MU) /(  V**2) 

RETURN 

END 

FUNCT  ION  PK(-'U) 

COMMON  AO , A  1 , A2 , A3 , 0 1 , B2, B3 , ALGi AL I N, ASO , ACURE , VO*PO , OELP , R3 , VI 0, 

+  V20,VP2,T02,T0D,TC1,T  2(1,  JC.PTT 
REAL  MU 

PK='1U*(8H-MU*(32*KU*«2)  ) 

RETURN 

END 

REAL  FUNCTION  LATliT  (T,P) 

LATHT  =  6 .26*4 . 18E-5 

RFTURN 

END 

/* 

A/GC .SYS  IN  DC  ♦ 

1.C34C-C4  11. M6  11.340  207.0  600. 0  293.0  96. 

2.7091  -2=9222  1.413  0.0  .54163  .749041  .605339 

0.005  1.0  3 

/* 

/* 
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IEF285I  SYSOUT 
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IEF285I  LISTER 

IEF285 I  VCL  SER  NOS  =  SCRQC1. 
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APPENDIX  B 


HUGONIOT  CURVE  OF  LIQUID  ARGON  OBTAINED 
BY  USING  THE  SIGNIFICANT  STRUCTURE 
MODEL  OF  LIQUIDS 

C.  T.  Tung 


I.  Introduction 

The  significant  structure  theory  of  liquids  has  been 
developed  by  Eyring  and  co-workers.  According  to  the  theory, 
a  liquid  is  considered  as  having  three  significant  structures; 
solid- like,  gas- like,  and  deganerate.  These  three  structures 
contribute  essentially  to  the  thermodyns.uic  properties  of  the 
bulk  system.  On  the  basis  of  these  considerations,  the  parti¬ 
tion  function,  f,  for  a  monatomic  liquid  such  as  argon  can  be 
* 

expressed  as'  . 


f  =  (a. 


3.3  /  83) 


N 


Vv 


(a|(V-V  )) 


N(V- 


vs)/v 


(N(V-VS)/V)! 


where: 

a1  =  exp(Eg/RT) 

&2  ~  1  -  exp(-8/T) 

83  =  1  +  ^h  exp(-s/RT) 

a4  =  (2rr>-rKT)^/h 


Superscripts  refer  to  literature  at  end  of  Appendix 


(1) 
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The  first  set  of  brackets  stands  for  the  solid-like  portion  of 
the  partition  function,  for  which  the  Einstein  oscillator  model 
is  used.  The  quantity  (1  +  n^  exp(-e/RT)  is  the  geometrical 
degeneracy  factor;  the  remaining  portion  is  the  gas- like  part. 
Using  Sterling's  approximation,  y!  «  (y/e)^  ,  Eq.  (1)  can  be 
rewritten  in  more  compact  form, 


3™s/V 


N(V-V  )/V 


(ala3/a2J)  5  <VeV/N)  S 


(2) 


In  Eq.  (2),  the  number  of  neighboring  positions,  n^  ,  is  equal 

to  n(V-V  )/V  ,  and  the  energy  needed  to  occupy  a  vacant  site, 

s  s 

s  ,  is  equal  to  aE„V  /(V-V  )  .  Both  n  and  a  are  propor- 

o  S  S 

tionality  factors,  E  is  the  energy  of  sublimation. 

At  high  pressures  and  temperatures  a  few  corrections 
are  necessary.  Hence,  Einstein  partition  function  in  Eq.  (2) 
should  be  replaced  by^ : 

a1(((l-g)/a2)  +  Sa4v£1/3)3 

where  g  *  exp(-A8/T)  ,  is  the  vibrational  quantum  number  and 
Vf  is  the  molar  free  volume  in  the  solid.  may  be  repre¬ 

sented  by: 

V£  =  (<Vs/N)2/3  "  (b/4N)1/3)3 

where  b  is  the  van  der  Waals  constant  and  b/4N  is  the  net 
molecular  volume.  With  these  corrections,  the  partition 
function  for  a  monatomic  liquid  at  high  pressures  and 
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temperatures  should  be  written  as: 

f  =  (a1((l-g)/a2)  +  ga5((Vs/N)1/3  -  (b/4N)1/3))3  (3) 

where : 

a5  =  1  +  b2b3 

b2  =  H(w-vs)/vs 

b3  =  exp(-b4) 

b4  =  aEsVs/RT(V-Vs)  . 

In  addition  to  the  corrections  mentioned  above,  the  pressure 

effect  on  V  must  also  be  considered.  For  moderately  high 
s 

pressures,  Eyring  suggests  the  linear  correction: 

V8‘  =  V8(l-Pp)  (4) 

In  order  to  extend  calculations  to  higher  pressures  we 
replace  (4)  by: 

vsf  =  vs  exp(-Pp)  ,  (5) 

which  reduces  to  Eq.  (4)  for  small  pressures  and  also  has 

positive  curvature,  which  is  necessary  for  shock  stability. 

Hence,  8  is  solid  compressibility  and  p  is  excess  pressure 

(2) 

above  a  standard  pressure. 

Krowing  the  totalpartition  function  as  a  function  of  T 
and  V  ,  we  are  able  to  calculate  thermodynamic  equations  of 
state  from: 
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A 

-  -kT  in  f 

(6) 

P 

-  -(5A/3V)T  - 

kT(5  in  f/av)T 

(7) 

E 

*  -T2(3(A/T)/3T)v  =  kT2(&  in  f/&T)v 

(8) 

S 

-  -(&a/st)v  = 

k  in  f  +  kT(3  in  f/ST)v 

(9) 

A, 

P,  E,  and  S  are 

the  Helmholtz  free  energy, 

pressure, 

internal  energy,  and  entropy. 

II.  Calculation  of  the  Hugoniot 
Curve  for  Liquid  Argon 

As  mentioned,  the  partition  function,  £,  is  a  function 
of  temperature  and  molar  volume.  Obviously,  to  calculate 
pressure,  internal  energy,  and  entropy  from  Eq.  (6)  to  Eq.  (9) 
is  straightforward,  but  tedious.  The  results  are: 

P  =»  Rt(l+B+C+((V-V  )/V2)  -  (3V  /V2)4n(a2F))  (10) 


E  =  RT  vD+J) 


(ID 


S  =  (E/T)  +  RTJ  +  (RVS/V)((ES/RT)  +  in  (a5/a23))  (12) 

where : 

L  =  -(Vs/V2)((Es/RT)  +  *n(a5/a23)) 

B  =  (b3VgV)(l  +  b4) 

c  =  (vs/v2)(h  +  in(T3/2V)) 

D  =  (Vs/V)((-Eg/RT2)  +  (30exp(-8/T)/T2a2) 

+  (b3TlaEs/a5RT2))  +  3(V-Vg)/2TV 
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F  =  ((l-g)/a2)  +  (vs1/3  -  (b/4)1/3)a4g/N1/3 

G  =  e(-Xg  +  l-a2  +  (A-l)  g  (l-a2))/(Ta2)2 

+  gT1/2((xe/T2)  +  (1/2T);  a4(Vs1//3  -  (b/4)1^3)  /  N1^3 

H  =  ^n(a43  e/N) 

J  =  (3Vg/V)((G/F)  -  e(l-a2)/T2a2) 

and  the  parametric  values  for  liquid  argon  are^3’^ 

n  (proportionality  factor)  =  10.7 

a  (proportionality  factor)  =  0.0052 

E  (sublimation  energy)  =  1888.6  cal/mole 

O 

b  (van  der  Waals  constant)  =  32.2  cc/mole 

!>  (vibrational  quantum  number)  =  5 

V  (molar  volume  of  solid  at  1  atm)  =  24.98  cc/mole 

o 

0  (Einstein  characteristic  temperature)  =  60.0  (°K) 

m  (atomic  weight)  =  39.944  gm/mole 

N  (Avogadro  number)  =  6.024  x  1023(mole)  ^ 

-27 

h  (Planck  constant)  =  6.6252  x  10  erg-sec 

0  =  2.5  x  10”3/atm  =  compressibility  of  solid 

argon  below  10^  atm. 

- 1  6 

k  (Boltzmann  constant)  =  1.380  x  10  erg/deg 

R  =  1.986  cal/raol  °K 

In  order  to  find  the  Hugoniot  cur^e  the  Rankine-Hugoniot  jump 
condition : 


E  -  Eq  =  2V?  +  P0)(V0  -  V) 


(13) 


is  required  to  be  satisfied.  The  problem  of  finding  the  Hugoniot 
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curve  is  equivalent  to  eliminating  both  T  and  E  among  Eqs.  (10), 
(11) ,  and  (13)  so  that  P  can  be  expressed  ir.  terms  of  V  only. 

The  relation  between  P  and  V  represents  the  Hugoniot  curve.  In 
principle  the  Hugoniot  curve  can  be  obtained  no  matter  how 
complicated  Eqs.  (10)  and  (11)  may  be.  But  in  practice  we 
accomplish  this  by  numerical  methods. 

The  calculation  procedure  is  schematically  as  follows: 

1.  Assign  a  value  for  V  (less  than  initial  volume  V  ) . 

2.  Guess  a  value  for  T  (higher  than  initial  temperature  T  ) . 

3.  Substitute  both  T  and  V  in  Ea.  (10)  and  calculate  the 
value  of  P. 

4.  substitute  the  calculated  P  into  Eq.  (5);  using  the 
new  V  (V  '),  P  is  recalculated  from  Eq.  (10).  This 

b  O 

process  is  repeated  until  consistency  is  obtained. 

5.  Use  present  T,  V,  and  V  to  calculate  E  from  Eq.  (11). 

O 

6.  Substitute  P,  V,  and  E  in  Eq.  (13)  which  can  be 
written  in  the  form, 

H(P,V,E)  =  E  -  Eq  -  *(P  +  Po)(Vc  -  V) 

If  H(P,V,E)  «  0,  P  is  the  right  value  which  corresponds 
to  the  assigned  V.  If  H(P,V,E)  «  0,  knowing  K  is 
positive  or  negative,  T  can  be  appropriately  adjusted, 
and  then  follow  with  step  (3).  This  loop  is  repeated 
until  H(P,V,E)  »  0  is  satisfied. 

From  this  double  iterative  method  the  relation  becween 
P  and  V  under  the  jump  condition,  Eq.  (13),  can  be  satisfied. 
This  relation  presents  the  Hugoniot  curve.  Furthermore,  using 
a  similar  method,  with  Eq,  (13)  replaced  by  Eq.  (12),  the 


6.5 

adiabatic  curve  can  also  be  obtained.  Details  of  this  double 
iterative  method  are  shown  in  the  computer  program  listing  at 
the  end  of  this  appendix. 

III.  Discussion  of  Results 

Each  isotherm  in  a  P-V  diagram  shows  both  the  existence 
of  a  maximum  pressure  and  a  discontinuity  at  V  *»  Vg  (See  Fig.  1) . 
This  is  due  to  using  the  Einstein  oscillator  model  for  the  solid 
partition  function.  In  the  Einsten  model  the  binding  energy  is 
assumed  to  be  volume-" independent . 

The  Huganiot  and  adiabatic  P-V  curves  can  exist  and 
have  been  calculated  only  in  the  region  well  to  the  right  of 
the  maximum-pressure  curve  of  Fig. 15  where  (dp/dv)^  <  0.  Both 
are  shown  in  Fig.  16  relative  to  the  isotherms.  In  Fig.  17  the 
calculated  curves  are  compared  with  measurements  reported  by 
van  Thiel  and  Alder.  The  agreement  is  remarkably  close  and 
suggests  that  minor  modifications  of  the  Eyring  theory  may  make 
it  valid  at  even  higher  pressures . 
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IV  G  LEVEL  0,  MOD  0 


MAIN 


DATE  =  6726* 


,VSC 


100 


500 

300 


FUGCMCT  CURVF  FOR  LIQUID  ARGON 
COMMCN/DATAM/M 
C0MMCN/ES02S/ES0 
COMMCN/APLCT/APi 101,  101) 

GQMMCN/CT VSE  T / T 1 ,T1 T,EXPT,DCWN,EXPST 
CUHMCN/YSET/Y ,OVCO 

-S^^n/>10sft/a’ESvVS,sni,ss*c,<*eeo»ppo»vo,to,ckk 
CQMMCN/SnC/i.LANK,  OOT,  C(  13) 

READ  CHARACTERS  IFOR  PLOTTING  CURVES) 

READ (5, ICC) BLANK, DOT, ( C{ K ) ,K* 1 , 1 3 ) 

FORMAT! 15A1) 

GIVE  PARAMETRIC  VALUES  FOR  LIQUID  ARGON 

A=0 .0052 

ES-1888.6 

£SO*ES 

VS-24.98 

SN=10.7 

SS=6C.O 

cK-I1I,-??i^^3;141*8*37,+n*5/,(6*625**3-0,’>,6*°2^.o)) 

LK- ALbG I CKK J  ♦  l  .0 

VSO=VS 

M.=*0 

READ  DOTH  INITIAL 
R£ AC ( 5  »3G0 ) V , T 


VOLUME  AND  INITIAL  TEMPERATURE 


$ 

*, 

c  IF  DATA  CARDS 

ARE  USEO  UP, STOP 

IF(V.LE.C.O) 

GO  TO  112 

£ 

M=NM 

£ 

WRITE(6,7C)M 

70  FORMAT! 'O',' 

DATA  MU 

SET  'CLANK'  IN  TWO  D! PENT IONAL  SPACE! 100*100) 

DO  210  K= 1,100 
DO  210  N-2,101 

ap«:k,n}-slank 

210  CONTINUE 

SST  UP  COORDINATES  8Y  'DOT' 

DO  211  K=  1 , 101 
AP  IK, 1 )  =  DCT 

211  CONTINUE 

DO  212  N= 1,101 
AP (1CI , N) =DOT 

212  CONTINUE 
TO=T 

vo*v 

CALCULATE  THE  ISOTHERMS 
CALL  PLOT  1 1  V, T ) 

CALCULATE  AN  A0IA8ATIC  CURVE  THROUGH  INITIAL  STATE 
CALL  ADIAB!V»T) 


lA/50/12 


IV  G  LEVEL  0,  MOO  0  MAIN  DATE  =  67264  14/50/12 

C  CALCULATE  THE  MJG0N1GT  CURVE 
CALL  SHCCK ( V  »T ) 

C  PLOT  THE  HUGCN1GT  CURVEiTHE  ADIABATIC  CURVEt  AND  THE  ISOTHERMS 
WRITE <6, 302)  ( (AP(IiJ) »J“1»1C1) »J-1»1C1) 

302  FORMAT  (  '  SICX.lOUl) 

VV= VC/39. 944 

C  WRITE  THE  INITIAL  STATE  UNDER  THE  HORIZONTAL  AXIS 
WRITE16.900) TO, VV 

900  F0RMAT(‘C*,3CX,' (INITIAL  STATE  T* • , F 7 . 2 t ' ( K )  V='.F10.5, 

1UCC/GM)  )•) 

C  GO  TC  RE AC  NEXT  DATA 
GO  TO  5C0 
112  CONTINUE 
RETURN 
END 
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G  LEVEL  0*  MOO  0 


DATE  *  67264 


14/50/12 


SUE  SHOCK 
FIND  HUGCN10T  CURVE 
SUBROUTINE  SHOCK!  V*  if} 

COMMCN/BCC/BLANKtOOT, C! 13) 

COMMON/ YSET/Y ,UVOO 

COMMON/C10SET/A,ES» VS*SN5SS»CK»EE0»PPCtVC*TC»CKK  , VSC 

VVaV/39.944 

WRITE (6»A5C) T ,VV 

450  FORMAT! •C,*2X»,HLG0N10T  CURVE • , • { I NI TI AL  STATE  T*'.F7.2. 
1  * l K )  V= ' tr 1C. 5f ' ICC/GM )  )•) 

■*RITE(6t201) 

201  FORMAT!  • C ' , 16X , • V ! CC/GM)  *,13X,'T!K)  M5X,  •  P  <B AR  )  • ,  12X t 
I'E  ICAL/GM) » , 13X» • SICAL/GM-K ) • , 10Xf »V/V0‘ ) 

VIN =V 
T  1N*T 

ASSIGN  A  CHARACTER  FOR  HUGONIOT  CURVE 
Y*C! 11 ) 

CALCULATE  ThE  INITIAL  VALUESIP*  AND  E» ) 

CALL  WRITE1IV.T) 

PPO=P ( V,T) 

EEO=EI V»T ) 

ASSIGN  A  VALUE  FOR  V  AND  FIND  THE  CORRESPONDED  P 
511  V=V-VO/25.0 
T  =  TO 
TL*T 

START  ON  THE  UOUeLE  ITERATIVE  CALCULATIONS 
IFIFl y,T).GE.O.O)  GO  TO  591 

SFTL=-l.O 
GO  TO  59? 

591  SFTL=l.O 

592  TU=TL+30C.0 
T*TU 

IF!F(V»T  ) . GE. 0 .0 )  GO  TO  593 

SFTU=-l.O 

GO  TO  594 

593  SFTU= 1.0 

594  IF($FTL*SFTU.LT.O.OJ  GO  TO  595 
TL=  TU 

SFTL=SFTU 

IF  ITU.GT. ICOOO.O)  GO  TO  611 
GO  TO  592 

595  TMMTL+TUJ/2.0 

UPPER  eOUNC  AND  LOWER  BOUND  HAVE  BEEN  FOUND 
T=TM 

IF ! F! V»  TJ.GE.O.O)  GO  TO  596 
SFTM=-1.C 
GO  TO 

596  $FTM=  1.0 


696 


V  G  LEVEL  Of  (■'CO  0 


SHOCK 


DATE  =  67264 


14/50/12 


696  IF(SFTN.EC.SFTU)  GO  TO  796 

TL=TF 

GO  TO  597 
796  TU=TP 

597  IF(  (TU-TU.GI.5.0)  GO  TO  595 
T= {TU+TL ) /2.G 

C  THE  SOLUTION  OF  T  IS  FOUNO 

C  USE  PRESENT  T  ANO  V  TC  CALCULATE  P 

CALL  WRITEl(VtT) 

1F( V.LT.VO/IO.O)  GO  TO  611 

GO  TO  511 
611  CONTINUE 
V= V  IN 
T=  T  IN 
RETURN 
END 


o  o  o 
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i  G  LEVEL  0,  HOD  0  MAIN  DATE  «  67264  14/50/12 


FIND  AN  ADIABATIC  CURVE  THROUGH  INITIAL  STATE 
(THE  METHOD  IS  SIMILAR  TO  THAT  FOR  FINDING  FUGCNIOT) 
SUBROUTINE  ADIAB(V.T) 

COMMON/ BCC/bLANK» DOT  *  C  { 13) 

C0MK0N/C10SET/A»ES» VS?SN»SS»CKjEEO»PPCt VO»TO,CKK  ,VSO 
CCMMCN/YSET/Y  »UVOO 
KRITE(6,450) 

450  FORMAT! 'C' t2X, 'ADIABATIC  CURVE') 
hR  IT  E( 6»  201 ) 

201  FORMAT  PC'  *loX,'V(CC/GM)  ',13X,'T(K)  • ,  15X ,  •  P  (BAR  )  • ,  12X, 
1' E (CAL/GN) ' » 1 3X, • S ( CAL/GM-K )'»10X»'V/VC') 

V  IN3  V 
T  IN=T 
Y=C ! 12 ) 

CALL  WRITEltVtT) 

S$0*S( V,T) 

511  V* V-V0/20  sO 
T=TO 
TL*T 

C  GET  THE  RIGHT  V* 

NEWVS=P ( V t  T i 

IF(S(V*T).GE. SSO)  GO  TO  591 

SFTL=- 1 .0 
GO  TO  592 

591  SFTL= 1*0 

592  TU=TLMOC.O 
T=TU 

NEKVS=P(V,T) 

I F ( S ( V *  T ) «GE<  SSO)  GO  TO  593 

SFTU='1.0 

GO  TO  594 

593  SFTU=l.O 

594  IF(SFTL*SFTU.LT.O.O)  GO  TO  595 
TL=TU 

sftl-sftu 

IF(TU.GT.ICOCO.O)  GO  TO  611 
GO  TO  592 

595  TM= ( TL+TU ) /2.0 
T  =  TM 

NEWVS=P(V»T) 

1  FIS (Vi T ) .GE.SSO)  GO  TO  596 

SF  TM=- 1 . 0 
GO  TO  696 

596  SFTM=*  l  .0 

696  IF(SFTM.EC.SFTU)  GO  TO  796 
TL*TN 
GO  TO 


597 
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G  LEVEL 

0,  HCO  0 

ADI  AS 

796 

TU*TH 

GO 

TO 

595 

597 

IF  (  (TU-TU.GT.5.0) 

T*( TL+TL ' /2. 0 

CALL  WRITEKV.T) 

IFIV.LT. V0/10. 01 

GO 

TO 

611 

611  CONTINUE 
V-V  IN 
T=T  I N 


RETURN 

ENO 


DATE  *  67264 


14/50/12 
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/  G  LEVEL  0,  POO  0  MAIN  DATE  *  67264  14/50/12 

C 

C  INTERNAL  ENERGY  FUNCTION 
FUNCTION  E(V,T) 

COHR0N/ClOSET/A,ES,VS*SN,SS,CK(E6O,PPO,VO»TO,CKK  ,VSC 

C0RRCN/CTVSET/T1 , T1 T t EXPT , DOWN . EXPST 

C0RRCN/E1FI/E1,F1,F12 

XL*5.0 

B=32.2 

E£VT*(VS/V)*<-ES/I2.C*T*m3.0*SS*EXPST/(T*T*(i.0-EXPSTm 
lSN*A*ES*EXPT/(2.0*T*T*D0WN))*1.5*I V-VS}/(T*V> 

E=2.C*T*T*EEVT 

F2l3($S/(T*T))<'(-XL*EXPST**XL  +  £X!’ST  +  {XL-l*0)*EXP$T**(XL  +  1.0)) 
l/l ( l.C-EXP$T)**2.0> 

F2aF2UF!2*!XL*SS/T*0.51/T 

El«(6„0*VS*T*7/V)*<F2/Fl-ISS*EXPST}/ (T*T*( 1 .O-EXPST  ) ) ) 

E*E*E1 

RETURN 

END 


o  o 
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«/  G  LEVEL 


0,  PCD  0  FAIN  OATE  =  6726A  14/50/12 

SUB  rvSET 

A  SET  OF  VARIABLES  IN  TERMS  OF  V  AND  T 
SUBROUTINE  TVSET(V.T) 

COKMOM/CIOSET/A.ES, VS, SN, SS.CK.EE O.PPC, VO, TC.CKK  ,VSG 
COFFCN/CTVSET/n.Tir,  EXPT.OCWN^EXPST 
Tl=A*ES*VS/(2.0*(V-VS>) 

T l T*T l/T 
EXPT=EXPI-T1T) 

00kN=1.0  +  SN*{V-VS}'»EXPT/VS 
EXPST=EXP(-SS/T) 

RETURN 

END 
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V  G  LEVEL  0,  XCO  0  MAIN  DATE  *  67264  14/50/12 

C 

C  PRESSURE  FUNCTION 
FUNCTION  P(V.T) 

COMMON/ ESCES/ ESQ 

C0HMCN/C1GSET/A,ES,VS,SN,SS»CK,EE0,PPC*V0.T0.CKK  ,VSO 

C0MMCN/CTVSET/T1,T1T,EXPT,D0WN,EXPST 

COMMON/ El  FI /El » FI , F12 

XL*5.0 

B0«32.2 

B£TA*2. 5/2421.0 

M*0 

N*0 

P0LC*0.0 
VS*V SO 

IF  l  V. LE. VS)  GO  TO  99 
ES«ESO 

C  OBTAIN  THE  RIGHT  VALUES  FOR  THOSE  VARIABLES  IN'CCMMON  CTVSET' 

CALL  TVScT(V.T) 

100  AA*-(  V$/(V*V)  )*<GS/(2.0*T}-3.0*AL0G{1.0~EXPSmALUG(CCWN) ) 
8B=(SN/V)*(  1.0*Tin*EXPT/OOWN 

CC*(VS/IV*V)  )*<CKU.5*ALCGm  *ALOG(V) J 
PPVT=AA+B8+CC*(V-VS)/< V*V> 

P»2.0*T*PPVT 

FI  1=1  l.0-EXPST**XL)/(1.0-EXPST) 
bo* VS/VSO 

F12=(EXPST**XL)  *  (CKK**0«  333)*(V$**0.333-{B/4.0)**0.333)*(T**0.5) 
F1=F1 1+-F12 

Pl=-(6.0*VS*T/(V*V))*AL0GIFl*I l.O-EXPST)) 

P=P+P1 

IF ( ABS {  P-POLC )  . L  T .  1  ,.0 )  GO  TO  901 

POLD-P 

M*M*1 

IF(M.GT.SO)  GO  TO  101 
C  START  ON  THE  ITERATIVE  CALCULATIONS 
VS* VSO*EXP( -bETA*P ) 

IF(V.LE.VS)  GO  TG  99 
90  CONTINUE 

ES=ESO*(VSO/VS)**Oj 333 
CALL  TVS  er  C  V , T  > 

GO  TO  100 
99  iN«N* l 
VS=V/2.0 

IFIN.GT.2)  GO  TO  101 
GO  TO  98 

101  CONTINUE 
901  CONTINUE 

RETURN 

END 


o  o  o 


81 


V  C  LEVEL  0,  MCD  0  FAIN  CATE  =  67264 


CALCULATE  THE  ISOTHERMS 
SLB  PLOT  1 

SUBROUTINE  PLOTKV.T) 

COMMON/OATAM/M 
COMMCfl/ 8CC/BL  ANK»  DOT  *  C  !  13) 

COMMON/ YStT /Y  tUVOO 
UR  I TE ( 6»  45C I 

450  FORMAT! 'C* ,2X,* ISOTHERMS* I 
V  IN»  V 
T  IN=  T 

IF1M.LE.21  GO  10  30C 

0  T=  50 . 0 
GO  TU  301 

300  0  T  = 1 0 0 . 0 

301  CONTINUE 
T=T-f;f 

00  1C  K=l,5 
C  ASSIGN  4  CHARACTER  FCR  PLOTTING 

Y=C!K  1 
T=T+CT 

wk i r e (b • a 5 1 )  k,t 

451  FORMAT! *C* ,2>  ISOTHERM  *, 1 2 , 5X ,* TEMPERATURE  • , F 1 0  .4 ,  •  t K  )  • 
wR  I TC- (  6 1  20 1 1 

201  FORMAT ( *C», lfcX, 'VICC/GM1 • , 1 3X , • T ( K) • , 1 5X , » P ( BAR  1 • , 12X , 

1  •  £  ( CAL/GM ) '  1 1 3X , '  S ! CAL/GM.-K I  •  ,10X,'V/V0*) 

V=  V  IN  *V  IN'/  10*0 
00  10  N=  1 ,  1 0 

V=V-V IN/IO.O 
C  CALCULATE  P 

CALL  V.R  I  T  E  l  <  V  »  T  l 
10  CONTINUE 
V=V  IN 

t  =  r  in 

RETURN 

END 


1 A/5C/12 
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V  G  LEVEL  0,  *00  0  MAIN  DATE  *  67264  14/50/12 

C 

SUBROUTINE  V.R1TE1 1 V*  T 1 
COMMON/ DAT AM/M 
COMMON/APLGI/AP ( 10 1 » 10 1 ) 

COMMCN/YSET/Y  »UVOQ 

COMMCN/CIOSET/A, ES» VS»SN»SSfCK*EEOt  PPO»VO* TO»CKK  , VSO 
VOO3 VO/39. 944 
T  IN*T 
V  IN*V 

C  CALCULATE  P,£,ANO  S 

PP3P(V, n 
EE*EC'v»T  J 
Sl*StV,T) 

C  CHANGE  UNITS 

V3V/39.944 

PP3PP*lOCG.O/( 24.21*0.987) 

EE=EE/39. 944 
Sl=Sl/39.944 
VR=V/VCO 

C  WRITE  ViT  tPtEfS  AND  V/V. 

WRITE(6,2C2)V,T,PP,EE.S1,VR 
202  FORMAT ( •  *,7X,6E20.6) 

C  LOCATE  THE  PCINT(P,V) 

UVOO=VCO/ICO.O 
N=1 .5+V/UVCO 

IFIM.EG.l)  GO  TO  300 

K=(  ICOOC.O-PP)  /100.0U.5 
GO  TO  in 
300  K=I2C.0-PP)/G.2*1.5 
111  CONTINUE 

C  GET  RIO  OF  ThOSE  POINTS  OUTSIDE  THE  REGI CN( 101*1  Cl ) 

IF (K.GT. 101 1  GO  TO  400 

IF ( K.LT. 1 )  GO  TO  400 
IF(N.GT.lOl)  GO  TO  400 

IF(N.LT.l)  GO  TO  400 
APIK,N)*Y 
400  CONTINUE 
V=VIN 
T=TIN 
RETURN 
ENO 


PART  B 


ACOUSTIC  WAVES  FOLLOWING  STRONG 
SHOCK  WAVES 
G.  R.  Fowles 

I.  Introduction 

Conventional  dynamic  equation  of  state  experiments,  in 
which  the  wave  and  particle  velocities  of  plane  shock  waves  are 
measured  in  a  sample,  yield  only  partial  information  about  the 
state  of  the  shocked  material.  This  information  comprises  the 
stress  component  normal  to  the  wave  front,  the  density,  and  the 
internal  energy.  In  particular,  the  normal  stresses  across  a 
plane  perpendicular  to  the  shock  front  are  not  determined.  Know¬ 
ledge  of  these  stress  components  in  addition  to  the  stress  normal 
to  the  front  is  tantamount  to  knowing  the  shear  modulus  and  the 
yield  strength  of  the  material  under  shock  conditions. 

In  an  elastic-plastic  solid  the  shear  modulus  and  yield 
strength  must  be  known  in  order  to  treat  problems  involving 
interactions  of  shock  and  rarefaction  waves;  a  simple  example  is 
that  of  a  decaying  shock.  Shock  attenuation  experiments  on 
aluminum  and  other  materials  have  shown  that  material  rigidity, 
characterized  by  the  yield  stress  and  shear  modulus,  has  a  sig¬ 
nificant  effect  on  shock  attenuation  at  pressures  up  to  at  least 

200  Kbar.  Moreover,  the  values  of  these  parameters  are  not 

(12  3") 

simply  predictable  from  known  zero-pressure  values.  5  ’  J 

Attempts  to  determine  the  shear  modulus  and  yield 
strength  by  means  of  one-dimensional  shock  attenuation  experi¬ 
ments  have  been  only  partially  successful.  Spallation  of  the 
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free  surface  on  which  measurements  are  made  severely  limits  the 
information  obtainable. 

In  this  paper  I  report  preliminary  results  of  a  study 
of  small  amplitude  wave  propagation  in  an  elastic-perfectly 
plastic  solid  considered  to  be  previously  stressed  to  the  yield 

point  in  uniaxial  strain . as,  for  example,  by  a  uniform  plane 

shock.  By  relaxing  the  restriction  that  the  flow  be  strictly 
one- dimens iona 1 ,  i.e.,  by  allowing  the  (plane)  acoustic  waves 
behind  the  shock  to  propagate  at  arbitrary  angles  with  respect 
to  the  direction  of  propagation  of  the  shock  one  finds  that  four 
distinct  acoustic  waves  are  possible,  compared  with  two  for  the 
one-dimensional  case.  Their  velocities  depend  in  general  on  the 
shear  modulus,  and  their  amplitudes  on  the  yield  strength.  Thus, 
there  are  a  greater  variety  of  measurements  possible  in  the  two- 
dimensional  case  than  in  the  one-dimensional  case.  This  result 
is  promising;  however,  it  is  not  yet  clear  how  best  to  make  use 
of  these  waves  experimentally.  They  can  be  generated  by  such 
means  as  reflection  at  interfaces  oriented  obliquely  to  the 
direction  of  shock  propagation. 

Another  application  of  the  theory  is  to  problems  such  as 
oblique  reflection  of  shocks  at  interfaces.  These  problems 
have  not  yet  been  investigated,  but  their  solution  is  a  natural 
extension  of  the  results  reported. 

II.  Fundamental  Relations  and 
Initial  Conditions 

As  the  starting  point  for  the  problem  we  assume  a  plane 
shock  propagating  in  the  x^  direction  in  an  isotropic  elastic- 
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perfectly  plastic  solid  satisfying  the  V.  Mises  yield  criterion. 
The  amplitude  of  the  shock  is  arbitrary  except  that  it  must  be 
large  enough  to  bring  the  material  to  the  yield  point.  Certain 
fundamental  relations  to  which  we  will  make  reference  are  listed 
below. 

We  assume  Cartesian  coordinates  x^(i  =  1,2,3)  and  let 
u^  be  the  velocity  of  the  material  at  point  x^.  Strains  are 
assumed  small  and  the  strain  rate  is  therefore  given  by: 


au. 

3 


For  elastic  strains,  Hooke’s  law  yields, 


U  *  V)  ay  -  v  akk  « 


(1) 


where  a^..  is  the  stress  tensor,  E  and  v  are  Young’s  modulus 
and  Poisson's  ratio,  and  e^  are  elastic  strains. 

The  plastic  strain  rate  tensor  is 


and  the  deviatoric  stress  tensor  is 


akk 


(2) 


The  second  invariant  of  the  deviatoric  stress  tensor  s,  is 
given  by, 


(3) 
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For  a  V.  Mises  solid ,  s  ^  k  where  k  is  the  yield  stress  in 
simple  shear — assumed  constant.  Plastic  deformation  occurs 
only  when 

s  =  k  ,  s  0  .  (4) 


Otherwise  the  deformation  is  elastic  and  n^j  =  e^.  .  When 
yielding  occurs  the  flow  rule  is 


t 


(5) 


where  K  is  an  undetermined  constant.  The  equation  of  motion 
(for  small  amplitude  waves)  is: 


a  ,  v 

aST  <°ij>  =  p  u  > 

where  p  is  assumed  constant. 

The  stress  matrix  of  the  original  state  is  diagonal, 
of  the  form: 


with  a22  =  a33  from  symmetry.  Moreover,  the  yield  criterion 
(Eq.  4),  written  in  terms  of  principal  stresses  is: 


(a^  -  o 2)  +  (<?2  ”  ^3)  ~  ^ 


6k2  =  2Y2 


where  Y  is  the  yield  stress  in  simple  tension. 
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For  matrix  A  this  implies 

~  a2  =  k  =  Y  . 

III.  Elastic  (Unloading)  Waves 
A.  Shear  Waves 

An  unloading  wave  by  definition  does  not  produce 

plastic  flow.  That  is,  s  <  0  .  For  these  waves,  Hooke's  law 

expresses  the  relation  between  stress  and  strain  rate  tensors, 

and  the  velocities  are  those  of  elastic  longitudinal  or  shear 

2  2 

waves,  i.e.,  Pc^  --  +  2\l  or  PC2  =  n  where  X.  and  n 

are  the  Lame  constants. 

We  assume  an  unloading  shear  wave  whose  wave  front  is 
parallel  to  the  axis  and  is  inclined  to  the  X-^  axis.  Its 
velocity  is  C2  =  «/n/p  ;  we  wish  to  find  its  amplitude  such 

that  superposition  of  the  stress  matrix  associated  with  the  wave 
and  the  initial  stress  matrix  just  maintains  the  material  at  the 
yield  point.  That  is,  we  wish  to  find  the  maximum  amplitude  of 
an  unloading  elastic  shear  wave  whose  wave  front  is  inclined  at 
angle  or  with  respect  to  the  original  shock  (Fig.  18). 

The  stress  matrix  associated  with  the  shear  wave  alone 
is 


~  bn 

b12 

0 

B  = 

b21 

b22 

0 

= 

bll 

b12 

0 

0 

0 

b21 

b22 

(8) 
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The  other  roots  are 
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Therefore  X  =  ±  J-Y  b^ 

The  eigenvectors  associated  with  the  shear  wave  can 
now  be  found  from 


' — 

— 

—  — 

—  — 

bll 

_b12 

b12 

'bu  - 

X1 

x2 

—  — J 

II 

X1 

x2  _ 

Expanding, 

bll  X1  +  b12x2  X  X1 

b12  X1  "  bllx2  X  x2  * 

Eliminating  •* 

x  X1  ■  bnxi  =  x  x2  +  bnx2 
x2  X1 


(10) 


Substituting  from  (10)  for  b-^ 


X  +  b 


11 

11 


*jnm 

Vx2;  1  +  X/Y 


(ID 


Equation  11  gives  the  principal  stresses  associated  with  the 
shear  waves  as  a  function  of  the  tangent  of  the  angle  between 
the  principal  axes  of  the  original  stress  matrix  and  the 
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principal  axes  of  the  shear  wave.  Figure  19  shows  a  plot  of 
this  relation.  The  ordinate  is  the  principal  stress  of  the 
shear  wave;  the  abscissa  is  the  angle  between  principal  axes 
(0) ,  or  the  angle  between  wave  fronts  (<*) .  (Note  that  the 
principal  axes  of  a  shear  wave  are  inclined  45°  with  respect 
to  the  wave  front.) 


B.  Longitudinal  Waves 

We  adopt  the  same  approach  as  for  shear  waves,  except 
that  the  matrix  B  is  now: 


B 


0 

0 


0 


b33 


Elastic  longitudinal  waves  are  characterized  by  the  relations 
between  principal  stresses: 


X2  X3 


1-v 


The  invariants  of  the  B  matrix  are  then: 


I1  =  bll  +  b22  +  b33  =  (t-t)  X1 

I2  =  "  (bllb22  +  b22b33  +  bllb33^+b122  = 
I3  =  bllb22b33  "  b33b122  =  (l^)2 


v(2-v) 

(1-v)2 


(12) 

(13) 

(14) 
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and  2  13  ^all+bll+a22+b22^  (16) 

±  ^/(a11+b11+a22+b22)2  +  4^bi22~^all+bll^a22+b22^* 
These  roots  must  satisfy  the  yield  criterion: 

!  1  o  I  I  9  l  I  o  9 

(x1-x2)z  +  (x2-x3) z  +  (x3-x1)z  =  2YZ  .  (17) 

Inserting  the  values  for  X^  from  (15)  and  (16)  into  (17);  we 
get  after  simplification: 

^all+bll^2  “  (aH+bH^  (a22+b22^  +  (a22+b22^2  (18) 

2  2  2 

*"  ( a 2^"b (3^ ^~tb l"^"a 2 2+b 22^  ( ^2 2"^"b 33)  3  b  12  ■  Y  . 


This  equation  is  to  be  solved  simultaneously  with  Eqs.  12,  13,  14 
to  establish  maximum  longitudinal  wave  amplitudes. 
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Recalling  that  an~a22  ~  Y  we  can  rewrlte  (18)  as: 

bll"4'b222+b332  +  Y(2bll_b22'b33^  "  bllb22  "  bllb33 


"b22b33  “3b12 


But  from  (13) 


blib22  +  bllb33+b22b33  b12  + 


o  v(l-v)  „ 
L  + - 7  X, 

d-v7  1 


Eliminating  b^  between  these  equations  and  employing  Eq.  12  to 
eliminate  b^  yields: 

1  r  ~  4-V *11  o 

3  bn  Y  “  xi  Y  +  1  - 9  !  X,  =  0  .  ( 

11  1-v  l  i.  (1_v)2  j  1 


We  now  note  that  is  a  principal  direction  of  both  matrices 

A  and  A'  =  A-*B .  (Eq .  15.)  Hence  b^  must  be  an  eigenvalue 
of  the  B  matrix,  i.e. , 

b33  =  X3  =  T^J  X1 
Hence,  from  (12)  again 


(20) 


Equations  19  and  20  give  b^  and  b22  in  terms  of  the 
principal  stress  of  the  longitudinal  wave,  parallel  to  the 
direction  of  propagation,  X^,  . 
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If  we  now  solve  for  the  eigenvectors  of  B  ,  we  get 


(bll'"X)xl  +  b12x2  " 


22  yA2 


(b33~X)x3  0 


x-,s2  X  -  b. 


0 


X  -  b. 


xi  ‘  r-v  +  bn 


i "  bn 


(X  =  Xn) 


and,  substituting  for  b^  from  (19), 


($2 


2v2  -  3v+l  -  (4v2-4v+1) 


4v  -  6v+2  +  (4v  -4v+l) 


This  expression  gives  the  amplitude  of  a  longitudinal  wave  whose 

-Ux  i\ 

wave  front  is  inclined  at  an  angle  tan  [ — j  with  respect  to 

x2 

the  original  shock  front  and  which  just  maintains  the  material  at 
the  yield  point.  It  is  of  the  form 


a  -  b 


2a  +  b 


where  a  and  b  are  functions  of  Poisson's  ratio,  v  . 

Figure  20  shows  a  plot  of  X^/Y  as  a  function  of  the  angle 


~  1  1 

a  =  cot  ( — )  for  several  values  of  v- „  Note  that  the  wave 

x2  1  o  X1  1 

is  compressive  for  angles  greater  than  about  55°  (—  =  . 

For  shallower  angles  the  wave  is  a  rarefaction  wave  and  can  be 

quite  large  for  large  values  of  Poisson's  ratio. 
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IV.  Plastic  (Loading)  Waves 

Loading  waves  occur  whenever  Eq.  4  is  satisfied.  The 
flow  rule  relating  the  strain  rate  tensor  to  the  stress  rate 
tensor  is  then  modified  by  the  addition  of  a  term  given  by 
Eq.  5.  Thus, 


v  • 

E  akk 


6 . .  +  Ko '  .  . 


We  take  the  direction  of  propagation  of  a  plastic  wave 
behind  the  shock  as  the  x^  direction  and  assume  the  other  axes 
are  oriented  so  that  the  shear  stresses  a.^;  and  cr^  ahead  of 
the  wave  are  zero.  That  is,  the  prime  set  of  axes  is  rotated 
about  Xn  with  respect  to  the  unprimed  set  (which  are  principal 
directions  o£  the  initially  stressed  material) .  The  wave  is  also 
assumed  to  have  infinitesimal  amplitude.  We  wish  to  find  the 
velocity  of  the  wave  as  a  function  of  the  angular  difference  be¬ 
tween  the  plastic  wave  and  the  initial  shock. 

Graggs  (5)  has  shown  that  an  infinitesimal  discontinuity 
in  stress  and  strain  propagates  under  these  conditions  with  a 
velocity  that  is  one  of  the  roots  of  the  quadratic 

ApCp4  -  Bp^Cp2  +  Cm-2  =  0  (21) 


where  Gp  is  the  plastic  wave  velocity,  n  the  shear  modulus,  P  the 
density,  and  A,  B,  and  C  are  given  by, 

A  =  (1  -  2v)k2 

B  =  (3  -  4v)k2  -  (1  -  2v)(o^2  -r  o122) 

C  =  2(1  -  v)k2  -  (1  -  2v)  c^2  -  2(1  -  v)au2 
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This  can  be  written  alternatively  in  terms  of  the  two  elastic 
velocities  and  C ^  and  the  bulk,  or  hydrodynamic  sound  speed, 
CH  ,  where 

r  2  _  X  +  2\x 

G1  ~  P 


K 

=  —  (K  =  incompressibility) 


In  terms  of  these  quantities,  Eq.  21  becomes, 


In  order  to  find  the  velocity  in  a  given  direction  we  need  to 

I 

know  the  stresses  ahead  of  the  wave,  anc*  °n  •  Taking 

a-^  and  &2  as  the  principal  stresses  behind  the  initial 
shock,  and  «  as  the  angle  of  rotation  of  the  prime  coordinate 
system  about  the  x,  axis  (i.e.  or  is  the  angle  between  the 
wave  fronts) ,  we  have 

all  =  al  cos  a  +  a2 

=  (-a^+a^)  sincr  cosa 


(2?) 
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Moreover , 


Hence , 


al  "  a2  Y 


2 

=  a-^  -  Y  sin  ct 


<7^2  =  “Y  sinct'coso 


p  -  Y ( 2/3  -  sin2<*) 


Inserting  these  values  into  Eq.  22  yields  two  velocities  for 
each  value  of  Poisson's  ratio  and  for  each  angle,  or  .  A  repre¬ 
sentative  case  is  plotted  in  Fig.  21. 

Several  features  of  these  curves  are  noteworthy.  For 
zero  angle  of  inclination  the  plastic  waves  travel  with  the 
velocities  of  hydrodynamic  and  elastic  shear  waves  respectively. 
For  other  angles  one  plastic  wave  speed  falls  between  C-^  and  CH, 
and  the  other  is  less  than  C2»  The  faster  wave  speed  increases 
to  at  just  that  angle  for  which  a  longitudinal  elastic  un¬ 
loading  wave  can  have  only  zero  amplitude  (Fig.  3).  These  waves 
are  in  general  mixed  waves  that  produce  changes  both  in  the  stress 
normal  to  the  front  and  in  the  shear  stress  tangential  to  the 
front.  Thus,  they  tend  to  rotate  the  principal  axes.  The 


stress  discontinuities  as  given  by  Craggs,  are: 

«  i  2 n(o  '  -I-  va  ') 


(a)  an  A  a. 


la33 


v  33  '  22  A 

(l-2v)  pc  2 


(23) 


98 


(b)  cr11  A  a , 


2i‘(°22  +  v  g33  l 
(l-2v)  PC  2 


}  A  a. 


(c)  A  a 


2  o 

®  1  o  /  c  C  -  \ 

7^- efr?) 4  ’u 

°11  vp  "C2 


It  is  easily  shown  that  for  any  plastic  wave  - -  >  0  ,  and  , 


moreover,  the  two  plastic  velocities  are  bounded  by  the  elastic 


velocities 


Co  <  cPl  <  c. 


Cp9  <  c. 


Hence,  from  the  equation  for  A  o ^  above  (23c),  we  see  that 
the  faster  wave  tends  to  decrease  the  shear  stress,  cr^  , 
while  the  slower  wave  tends  to  increase  it. 

Some  appreciation  for  the  structure  of  a  finite  amplitude 
plastic  wave  can  be  gained  by  numerical  integration  of  the  above 
equations.  Each  infinitesimal  wave  front  alters  the  stress 
state  behind  it  and  since  the  wave  velocity  depends  on  the 
stress  state  ahead,  finite  amplitude  waves  will  generally  show 
amplitude  dispersion. 

Fig.  22  shows  a  plot  of  the  stress  normal  to  the  (fast) 
plastic  wave  as  a  function  of  the  wave  velocity  for  particular 
values  of  Poisson's  ratio  and  the  angular  difference  between  the 
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plastic  wave  and  the  initial  shock.  The  angular  difference  in 
this  case  is  such  that  Che  head  of  the  plastic  wave  travels  with 
velocity  C-^.  In  general  it  is  slower  than  so  that  there  is 
a  region  of  uniform  stress  between  elastic  and  plastic  wave 
fronts.  As  the  wave  speed  approaches  hydrodynamic  wave  speed, 
increasingly  large  increments  in  the  stress  normal  to  the  front 
are  required  for  a  given  increment  in  wave  speed.  Thus,  the 
shear  stress  only  asymptotically  tends  to  zero  and  the  wave  speed 
asymptotically  approaches  hydrodynamic  speed. 

This  model  does  not  permit  the  formation  of  a  shock  front 
as  a  true  discontinuity  in  stress  although  the  stress  gradient 
becomes  larger  with  increasing  stress.  However,  we  recall  that 
the  equations  are  based  on  the  assumption  of  small  amplitudes 
and  the  equation  of  motion  therefore  has  no  convective  term.  It 
would  be  of  interest  to  extend  the  theory  to  include  finite 
amplitudes . 
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FIGURE  CAPTIONS 

Fig.  IB.  Wave  front  configuration. 

Fig.  19.  Maximum  shear  wave  amplitudes  as  function  of  angle 
of  inclination  of  wave  fronts. 

Fig.  23.  Maximum  diiatational  wave  amplitudes  as  function 
of  angle  of  inclination  of  wave  fronts. 

Fig.  21.  Plastic  wave  velocities  as  function  o^  angle  of 
inclination  of  wave  fronts.  Poisson’s  ratio, 
v  =  0.25. 

Fig.  22.  Normal  stress,  ,  of  plastic  wave  as  function  of 
velocity.  Poisson's  ratio,  v  =  0.30,  inclination  of 
wave  fronts,  «  =  51.6°  =  0.9  radian. 
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Shear  Wave  Amplitudes  as  Function  of  Angie  of  Inclination  of  Wave 
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Fig.  21 

Plastic  Wave  Velocities  as  Function  of  Angle  of 
Inclination  of  Wave  Fronts 
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CP/CH 
Fig.  22 

Normal  Stress,  <y,,  ,  of  Plastic  Wave  as  Function  of  Velocity 


PART  C 


PHASE  TRANSITIONS  UNDER  DYNAMIC  CONDITIONS 
M.  H.  Miles 

I.  Physical  Considerations 

The  purpose  of  this  paper  is  to  survey  our  present 
understanding  of  shock  induced  phase  transformations  in  solids. 
It  is  commonly  observed  that  a  stable  crystal  structure  at  a 
given  temperature  and  pressure  becomes  unstable  upon  change  of 
temperature  or  pressure.  In  general  we  expect  an  increase  in 
pressure  to  favor  rearrangement  of  the  atoms  into  a  crystal 
structure  that  minimizes  the  volume  while  an  increase  in  tem¬ 
perature  Favors  an  arrangement  of  atoms  that  maximizes  entropy. 

From  the  most  fundamental  viewpoint  we  would  desire 
to  be  able  to  predict  in  advance  the  equilibrium  structure  crom 
the  known  properties  of  isolated  atoms.  This  x-?ould  entail  the 
many-bodied  quantum  mechanical  calculations  of  the  cohesive 
energy  for  the  various  likely  crystal  structures  as  functions 
of  pressure  and  temperature,  unfortunately  the  cohesive  energy 
for  the  different  crystal  structures  even  far  removed  in 
temperature  or  pressure  from  a  phase  boundary  are  not  too 
different.  The  difference  is  often  of  the  order  of  a  percent 
of  the  cohesive  energy  and  this  is  within  the  accuracy  of  the 
quantum  mechanical  calculation.  It  is  doubtful  if  the  fundamen¬ 
tal  approach  is  capable  of  reliably  predicting  in  advance  the 

phase  boundaries  on  a  PVT  diagram. 
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Another  less  fundamental  approach  is  to  avoid  the 
quantum  mechanics  by  using  a  model  calculation  based  on  an 
assumed  atomic  interaction.  Again  the  details  of  the  calcula¬ 
tion  are  complicated  and  only  approximate  solutions  to  the 
model  are  obtainable.  Even  if  the  model  seems  to  favor  a 
given  crystal  structure  it  is  always  somewhat  uncertain  if  this 
is  due  to  the  mathematical  approximations  or  to  the  model  itself. 
Since  effective  forces  between  atoms  extend  beyond  nearest 
neighbor  as  evidenced  by  differences  in  cohesive  energy  between 
hexagonal  close  packed  and  face  centered  cubic  structures,  the 
more  tractable  models  are  expected  to  be  poor  approximations  to 
the  real  binding  forces. 

A  third  approach  is  an  engineering  one  based  upon  some 
thermodynamical  model  for  the  phase  transformation.  Often  this 
approach  is  a  combination  of  empirical  and  semi-empirical  cor¬ 
relations  largely  based  upon  experimental  observations  of  a 
particular  system.  A  phase  transition  is  feasible  if  the  Gibbs 
free  energy  for  the  rival  structures  are  equal.  We  may  indicate 
this  situation  by  showing  the  Gibbs  free  energy  and  enthalpy  for 
two  polymorphs  as  functions  of  temperature  at  constant  pressure 
and  by  showing  the  Gibbs  and  Helmholtz  free  energies  of  the  two 
polymorphs  as  functions  of  pressure  at  constant  temperature. 
Confidence  in  this  approach  necessitates  accurate  calculations  of 
the  required  thermodynamic  quantities  or  at  least  a  useful 
representation  of  the  necessary  thermodynamic  quantities 
empirically  from  experimental  observations.  A  partial  solution 
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or  at  least  a  useful  classification  of  phase  transformation  is 
obtained  with  the  aid  of  the  Clapeyron  equation.  The  Clausius- 
Clapeyron  equation  gives  the  variation  of  the  transition  pressure 
with  temperature.  This  requires  accurate  calculation  of  the 
volume  change  and  entropy  change  but  useful  information  is 
obtained  if  merely  the  sign  of  the  entropy  and  volume  changes 
can  be  obtained. 

In  attempting  to  understand  shock-induced  phase  trans¬ 
formations  one  would  like  answers  to  such  basic  questions  as 
(1)  why  does  the  transformation  occur,  (2)  what  is  the  mechanism 
for  the  transformation,  and  (3)  what  differences  exist  in  static 
versus  shock  induced  transformations.  We  have  touched  upon  some 
of  the  difficulties  pertinent  to  the  first  question.  In  dynamic 
shock  the  shifting  of  atoms  to  a  new  structure  must  occur  in 
times  of  the  order  of  the  transient  pressure  duration.  This  time 
is  short,  being  of  the  order  of  microseconds.  The  situation 
under  dynamic  shock  may  well  be  influenced  more  by  the  kinetics 
than  by  the  equilibrium  thermodynamics  of  the  transformation. 

Since  our  present  fundamental  understanding  of  the 
stability  of  solids  is  in  a  rudimentary  state,  it  seems  advisable 
to  carefully  review  the  present  experimental  knowledge  concerning 
pressure  induced  phase  transformations.  The  "130-kbar"  poly¬ 
morphic  transition  in  iron  has  attracted  considerable  static  and 
dynamic  experimentation.  The  first  observation  of  this  polymorph 
was  reported  in  1956  in  the  shock-wave  investigation  of  Bancroft, 
Peterson,  and  Minshall.^  For  some  time  the  nature  of  the 
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transformation  was  considered  to  be  the  bcc  to  fee  (a  to  y) 
transformation  of  iron  such  as  occurs  at  910°C  and  atmospheric 

rt  n  / 

pressure.  Several  investigators^ *  *  investigated  the  «  to  y 
phase  boundary  out  to  pressures  of  about  90  kbar.  The  tempera¬ 
ture  dependence  of  the  phase  line  around  room  temperature  was 
investigated  by  Minshall'5  under  dynamic  conditions  which  seemed 
inconsistent  with  any  reasonable  extension  of  the  ot  to  y  phase 
line.  This  led  Fowler,  Zukas,  and  Mi.ishall^  to  question  the 
alpha  to  gamma  transition  suppositior. .  Johnson,  Stein,  and 
Davis^  reported  in  1962  shock  compression  results  on  specimens 
in  the  temperature  range  of  70°K  to  1158°K.  For  temperatures 
up  to  about  500°C  results  similar  to  Bancroft's  were  obtained 
while  above  500°C  a  transition  that  was  much  more  pressure 
dependent  was  indicated  in  fair  agreement  with  the  low  pressure 
ff  to  y  statically  determined  phase  boundary.  Johnson's  et  al. 
temperature-pressure  data  together  with  microstructual  observa¬ 
tions  suggested  a  triple  point  at  about  110  kbar  and  500°C. 

They  concluded  that  the  low  pressure,  high  temperature  phase  line 
was  the  alpha  to  gamma  transition  while  below  500°C  the  trans¬ 
formation  was  from  ot  to  an  "x"  phase  different  from  y,  being 
most  likely  hep. 

O 

Balchan  and  Drickamer  obtained  the  phase  change  stat¬ 
ically,  observing  a  sharp  change  of  resistance  at  133  kb  and  20°C. 
The  first  high  pressure  X-ray  investigation  was  performed  by 

Q 

Jamieson  and  Lawson.  In  the  high  pressure  phase  region  at  room 
temperature  they  observed  an  extra  X-ray  line  that  agreed  with 
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known  volumes  significantly  better  if  it  were  assigned  to  an 

intense  hep  line  rather  than  to  a  corresponding  fee  line. 

Later  improved  X-ray  work  by  Clendenen  and  Drickamer^  and  by 

Takahashi  and  Bassett^  established  that  the  high-pressure 

phase  at  roam  temperature  was  indeed  a  hexagonal  structure. 

12 

More  recently  Bundy  has  confirmed  the  pressure- temperature 
phase  diagram  for  iron  placing  the  triple  point  at  110  ±  3  kbar 
and  490  ±  10°C.  Bundy  calibrated  his  data  by  assuming  that  the 
110  kbar  assignment  by  Johnson  et  al.  for  the  triple  point  was 
correct.  The  justification  being  that  their  room  temperature 
data  correlated  well  with  the  130  kbar  a,  e  shock  transition 
of  Bancroft. 

13 

Loree  et  al.  has  recently  studied  the  dynamic  trans¬ 
formation  for  pure  iron  obtaining  for  the  best  value  for  the 
onset  of  the  dynamic  transition  129  ±  1  kbar  at  room  temperature. 
Earlier  transformation  pressures  appear  to  be  too  high  for  two 
reasons:  (1)  the  possibility  of  overdriving  the  transformation 

with  excessively  high  input  pressures,  and  (2)  the  samples  were 

/ 

not  annealed.  It  is  expected  that  annealed  samples  would  have 

the  lox^est  elastic  wave  and  therefore  the  lowest  transition 

pressure.  This  seems  to  be  indicated  by  the  work  of  Loree  et  al. 

1  /. 

It  is  interesting  that  Bundy"  has  presented  evidence  that 
dynamic  and  static  pressures  for  initiation  of  transformations 
are  identical  for  pure  iron  but  the  static  pressures  for  iron 
alloys  of  V  and  Co  show  much  larger  increases  compared  to  the 
dynamic  pressures  as  the  percentage  of  V  or  Co  is  increased. 

The  difference  at  20  wt  %  Co  is  huge > being  about  288  kb  statically 
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compared  to  136  kb  for  shock. 

Another  transformation  with  considerable  static  and 
dynamic  experimentation  occurs  in  bismuth.  Duff  and  Minshall^ 
were  the  first  to  observe  shock  induced  phase  change  in  bismuth. 
Their  shock  data  for  specimen  temperatures  of  -27°,  42°,  87°, 
and  236°C  indicated  a  transition  about  3.5  kbars  higher  than 
the  statically  determined  phase  diagram  of  Bridgman.  The  slope 
of  the  shock  data  was  -50.8  bars/°C  compared  to  the  statically 
determined  slope  for  the  Bil  to  Bi  II  phase  line  of  -50  bars/°C. 

On  this  basis  it  was  assumed  that  the  high  pressure  phase  was 
Bi  II  even  though  the  samples  were  subjected  to  shock  pressures 
far  into  the  Bi  III  static  equilibrium  region.  The  236°C  shocked 
crystal  was  driven  into  the  liquid  bismuth  region  of  the 
equilibrium  phase  diagram.  Since  melting  is  considered  to  be 
a  slow  process  compared  to  shock  pressure  durations  and  there 
apparently  was  no  evidence  for  melting,  it  appears  that  this  is 
a  clear  example  of  a  shock- induced  transition  to  a  thermodynamically 
unstable  crystal  lattice  instead  of  to  the  stable  liquid  phase. 
Larson‘S  has  repeated  room  temperature  shock  investigation  of 
bismuth.  Larson  measured  the  sample  pressure  using  quartz 
pressure  gauges  whose  readings  were  calibrated  assuming  linear 
Hugoniots  for  bismuth  and  quartz.  After  adjusting  the  observed 
dynamic  transition  pressure  to  an  effective  hydrostatic  pressure 
Larson  achieved  a  transition  pressure  of  25.4  kbars  for  isotropic 
bismuth  and  25.9  kbars  for  large  grain  cast  bismuth.  There  was 
no  overdriving  of  the  transition  pressure  e^en  for  samples  down 
to  1.5  mm  in  thickness.  Since  the  transit  time  of  the  shock 
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wave  through  such  a  thin  sample  is  less  than  a  microsecond,  the 

characteristic  time  for  the  transformation  is  much  less  than  a 

microsecond, being  perhaps  of  the  order  of  a  few  nanoseconds. 

The  static  transformation  pressure  at  room  temperature  has  been 

determined  to  excellent  precision  by  Kennedy  and  La  Mori^  to 

be  25.4  ±  0.1  kbar.  It  appears  that  bismuth  and,  perhaps  also 

pure  iron, cannot  be  overdriven  even  for  very  thin  specimens. 

The  importance  of  sizeable  shear  stress  in  reducing  the 

nucleation  and  growth  times  is  suggested  by  comparing  the  shock 

results  with  the  pure  hydrostatic  pressure  results  of  Davidson 
18 

and  Lee.  Delay  times  for  initiation  of  the  high  pressure 
phase  of  the  order  of  several  minutes  were  observed  for  both 
poly  and  single  crystal  bismuth  followed  by  slow  growth  of  the 
high  pressure  phase.  The  transition  pressure  and  transforma¬ 
tion  rate  were  found  to  be  independent  of  the  presence  of  grain 
boundaries.  It  seems  that  for  very  low  shear  stresses  and 
pressures  only  slightly  above  the  transition  pressure  that  the 
transformation  favors  thermally  activated  nucleation  and  growth 
processes.  However  in  the  shock  data  for  bismuth  the  new  phase 
must  nucleate  extremely  fast  and  the  new  phase  must  propagate 
in  the  shocked  sample  with  a  velocity  close  to  the  sound  velocity. 
Other  materials  such  as  antimony  relax  into  a  new 

structure  much  more  slowly  than  iron  or  bismuth.  Minshall's 

19 

work  on  antimony  referred  to  by  McQueen  showed  overdriving 

of  the  transition  pressure  for  samples  thicker  than  20  mm. 

20 

Warnes  has  recently  confirmed  and  extended  Minshall's  earlier 
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work.  Apparently  the  overdriving  is  due  to  delay  in  nucleation 
or  an  initial  slow  growth  process.  It  is  reported  by  the  data 
of  Breed  and  Venable‘S  from  the  PHERMEX  facility  that  X-ray 
photographs  show  that  the  plastic -two  wave  forms  at  the  sample 
interface  and  accelerates  rather  slowly  to  its  characteristic 
velocity.  This  gives  a  time  dependent  phase  transition  with 
the  plastic  2  wave  being  delayed  about  0.6  microseconds. 

It  is  obvious  that  the  short  duration  of  the  transient 
pressure  pulse  places  severe  limitation  on  any  mechanism  of 
transformation  that  requires  appreciable  time.  This  suggests 
that  shock- induced  transformation  should  be  considered  to  be 
classified  as  Martensitic  among  the  vast  literature  of  solid 
transformations.  Certainly  any  growth  by  diffusion  of  atom  by 
atom  across  the  interface  simply  requires  orders  of  magnitude 
too  much  time.  The  individual  atoms  must  undergo  a  correlated 
relative  movement  of  somewhat  less  than  one  interatomic  distance. 
This  correlated  atomic  shuffles  or  movements  are  similar  to  what 
occurs  for  example  during  mechanical  turning.  There  are  many 
examples  of  temperature  induced  Martensitic  transformations 
that  are  fast  enough  to  suggest  that  similar  atomic  shuffles  are 
initiated  by  pressure  pulses.  It  is  felt  that  a  study  of  the 
Martensitic  transformations  will  shed  light  upon  the  transforma¬ 
tion  process  and  that  shock  studies  may  well  prove  a  useful 
approach  in  understanding  the  martensitic  transformations. 

The  most  obvious  characteristic  of  the  martensitic 
transformation  is  the  so-called  shape,  deformation.  This  reveals 
itself  in  rather  well-defined  surface  distortions.  These  surface 
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relief  effects  usually  indicate  that  straight  lines  in  the 
crystal  are  transformed  into  straight  lines  and  planes  are 
transformed  into  planes.  It  is  also  known  that  the  martensitic 
phase  even  though  it  has  a  different  crystal  structure  has  a 
definite  lattice  orientation  relationship  to  the  parent  phase. 

The  particular  plane  of  the  parent  structure  called  the  habit 
plane  separates  the  two  phases.  For  convenience  we  will  make 
two  classifications  of  martensite  transformations.  The  most 
common  is  perhaps  the  platelike  martensite  which  forms  from 
numerous  nuclei  in  a  crystal  with  each  plate  apparently  growing 
independently  into  a  distinct  plate.  There  is  also  a  "single¬ 
interface  type"  martensite  which  occurs  in  some  materials  such 
as  Au-Cd  alloys.  In  a  single  crystal  the  parent-product  inter¬ 
face  extends  completely  across  the  crystal  so  that  the  interface 
plane  does  not  experience  the  volume  constraints  present  for 
platelike  martensite  formation.  The  boundary  between  the  parent 
phase  and  the  region  of  product  phase  is  planar  for  a  single 
crystal.  Included  regions  of  platelike  martensite  are  usually 
lenticular  in  shape.  The  shape  deformation  of  martensitic 
plates  constrained  by  the  parent  matrix  gives  rise  to  strain 
energy  that  may  be  very  large  so  that  further  growth  is  stopped. 
Additional  growth  upon  cooling  does  not  begin  until  the  chemical 
driving  force  can  overcome  this  strain  energy.  There  may  be 
competing  nucleaticn  and  growth  processes  which  begin  at  smaller 
driving  forces  giving  rise  to  the  oft  observed  martensite 
appearing  only  during  rapid  cooling  from  above  the  transformation 
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temperature.  The  reaction  starts  at  a  characteristic  temperature 
(Mg)  which  depends  upon  previous  mechanical  and  thermal  history 
and  on  grain  size.  For  martensite  in  steels,  the  chemical  driving 
force  is  about  300  cal/mole  but  for  other  solids  with  smaller 
shape  change  the  driving  force  may  be  smaller.  In  general  a 
large  driving  force  implies  a  large  temperature  hysteresis  between 
Mg  for  the  cooling  transformation  and  for  the  reverse  trans¬ 
formation  upon  heating. 

If  the  chemical  driving  force  is  not  large  enough  for 
spontaneous  transformation  or  even  of  the  wrong  sign,  martensite 
may  sometimes  be  produced  by  externally  applied  stress.  The 
lattice  transformation  may  be  viewed  as  a  mode  of  mechanical 
deformation  comparable  with  mechanical  twinning.  The.  shape  of 
the  mechanical  twins  are  often  very  similar  to  the  shape  of  mar¬ 
tensitic  plates. 

The  crystallographic  theory  of  martensitic  transformations 

no 

as  developed  by  Wechsler,  Lieberman,  and  Read'  and  a  fundamentally 

23 

equivalent  theory  by  Bowles  and  Mackenzie  is  essentially  phe¬ 
nomenological,  concerned  only  with  the  crystallographic  features. 
The  problem  of  nucleation  and  kinetics  remains  essentially  un¬ 
solved.  The  central  aspect  of  the  crystallographic  theory  is  to 
describe  the  proper  initial  and  final  atom  positions  and  to  satisfy 
experimentally  observed  shape  deformations  with  an  undistorted, 
unrotated  habit  plane.  Using  methods  of  matrix  algebra  it  is 
possible  to  transform  one  crystal  structure  into  another  but  an 
additional  matrix  is  generally  needed  to  give  the  correct  shape 
deformation  and  the  invariant  habit  plane.  P  lysically  the 
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additional  matrix  represents  twinning  or  dislocation  slip.  As 
far  as  the  theory  is  concerned  there  is  no  preference  given  to 
the  order  of  events  between  the  lattice  deformation  and  the 
crystal  deformation  or  between  the  choice  of  slip  or  twinning. 

Under  shock  conditions  it  is  not  known  if  the  phenom¬ 
enological  theory  applies.  Both  elements  of  a  crystal  deforma¬ 
tion  and  lattice  deformation  are  expected  to  exist  under  shock 
conditions  since  the  solid  has  been  driven  into  a  region  of 
plastic  relaxation  prior  to  relaxation  into  the  high  pressure 
crystal  structure.  The  restrictions  of  a  specific  shape 
deformation  and  an  undistorted,  unrotated  habit  plane  may  or 
may  not  remain  for  shock  conditions. 

In  discussing  the  kinetics  of  martensitic  transformations 
it  is  often  the  nucleation  rather  than  growth  that  is  rate 
determining.  The  work  of  Bunshah  and  Mehl^  on*  Jron-nickel- 
carbcn  alloy  indicates  that  the  linear  growth  of  individual 
plates  is  about  one- third  the  velocity  of  sound  in  the  alloy. 

The  velocity  was  observed  to  be  independent  of  temperature  in 
the  range  -20°  to  -200°C  indicating  that  the  growth  was  not 
thermally  activated.  This  interpretation  does  help  in  under¬ 
standing  athermal  martensite  where  the  nucleation  rate  is  a 
function  of  temperature  independent  of  time  and  the  understanding 
of  isothermal  martensite  where  the  nucleation  rate  for  a  partic¬ 
ular  temperature  is  time  dependent.  In  shock,  if  there  is  a 
delay  in  nucleation  of  the  high  pressure  phase  we  should  expect 
a  high  pressure  precursor  to  the  relaxation,  due  to  phase  change, 
which  decays  at  a  rate  depending  upon  the  nucleation  delay  time. 
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Warnes  has  suggested  that  as  some  of  the  material  behind  the 
shock  begins  to  relax  into  its  higher-density  form,  rarefaction 
waves  are  emitted.  The  forward  rarefaction  overtakes  the  over¬ 
driven  shock  thus  attenuating  it.  The  progress  of  the  plastic  2 
wave  is  being  delayed  by  the  relaxing  material  ahead  of  it. 

When  the  nucleated  region  relaxes  to  a  state  on  the  Hugoniot 
near  the  transition  pressure, emission  of  further  rarefaction 
is  no  longer  possible.  The  plastic  2  wave  is  presumed  to  now 
proceed  with  its  characteristic  velocity. 

At  present  we  have  little  or  no  understanding  why  a 
given  phase  transformation  behaves  as  observed.  What  seems 
totally  lacking  is  any  detailed  plausible  models  for  initiation 
of  a  new  phase  and  the  subsequent  kinetics  of  growth  yielding 
the  plastic  2  wavefront.  The  role  of  crystal  defects  in 
nucleation  of  the  new  phase  seems  so  far  essentially  unexplored. 
There  is  need  of  further  data  especially  on  the  simpler 
solid  state  systems.  For  example  the  simplest  martensitic 
transformation  is  from  a  high  temperature  fee  phase  to  a  low 
temperature  hep  phase  as  found  in  cobalt.  An  evaluation  of 
plausible  mechanisms  for  phase  transformations  seems  in  order. 
Useful  theoretical  proposals  should  be  amenable  to  experimental 
evaluation  and  sufficiently  realistic  to  be  taken  seriously  by 
shock  wave  experimenters. 

II.  A  Simple  Martensitic  Model 

Suppose  that  each  grain  of  mean  diameter  d  has  NQ 
nucleation  sites  distributed  around  its  boundary  and  that  N(P) 
of  these  are  activated  at  pressures  less  than  P.  Suppose 
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further  that  once  a  site  is  activated,  it  generates  a  platelet 
of  voli.me  <*d  which  runs  across  the  grain  in  time  a/c,  where 
c  is  sound  velocity  in  the  original  material,  Fig.  23. 

For  simplicity  take  the  densities  of  old  and  new  phases 
to  be  the  same.  Then  when  the  first  platelet  runs  across  the 
grain  it  transforms  a  mass  fraction  ad/d  to  the  new  phase. 

As  transformation  proceeds,  the  amount  of  mass  transformed  by 
activation  of  each  new  site  is  reduced.  Assume  that  when  a 
fraction  X  has  been  transformed,  activation  of  a  new  site 
increases  X  by  an  amount  (1-X)  ff/d  .  Then  for  very  slow 
increase  in  pressure,  X  can  be  assumed  to  equal  its  equilibrium 
value,  Xeq.  Then 

c 

dXeq  dN  (1-Xeq)<* 

__  = 

or 

Xeq  =  1  -  expC-oN/d^)  (2.1) 

A  graph  of  Xeq  vs  P  might  have  the  general  features  shown  in 
Fig.  24.  At  r? oiks  pressure  P-^  one  would  say  the  transforma¬ 
tion  started.  At  ?2  ^  would  be  effectively  completed. 

If  we  now  forego  the  earlier  assumption  that  the  two 
phases  have  equal  densities,  the  curve  of  Fig. 24  can  be  converted 
to  a  P-V  curve,  as  in  Fig.  25.  Here  OAC  and  QBD  are  compression 
curves  of  phases  1  and  '2,  respectively.  OAB  is  the  curve  obtained 
by  plotting  the  average  specific  volume,  v,  against  P,  where 

83  (1  -  X)v^  +  IV2 


v 


(2.2) 


*  - 
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and  X  is  taken  to  be  Xeq  of  Eq.  (2.1)  and  Fig.  24  The  slope 
of  OAB  at  any  point  can  be  determined  by  differentiating 
Eq.  (2.2)  and  combining  with  Eq.  (2.1): 


dv  dv,  dv«  dN 

—  -  (1  -  *  )—  +  *  —  +  (v2"vP  “7  —  ~  *  ) 

dP  eq  dP  eq  dP  1  1  dz  dP  eq 


(2.3) 


If  dv^/dP  »  dv£/dP  and  -  v-j  =  Dv  =  constant,  this 
reduces  to 


dv/dP  =  (dVl/dP)  +  Dv  »(1  -  Xeq)(dN/dP)/d2 


(2.4) 


In  its  integrated  form,  Eq.  (2.2)  is 


Vo  +  (v,  -  Vo)  e 


-tfN/d' 


(2.5) 


If  pressure  increases  very  rapidly  the  growth  of  plate¬ 
lets  may  fall  behind  the  pressure  increase,  so  X  may  have 
other  than  its  equilibrium  value .  Suppose ,  in  Fig  25 ,  that  the 
system  is  at  some  point  A  and  that  P  is  suddenly  increased 
to  P  +  6P  at  B.  Then  the  number  of  activated  sites  is  in¬ 
creased  to  N  +  (dN/dP)6p  and  X  starts  to  increase  toward 
point  C  at  the  rate 


dX 

dt 


f  <^>3 


(X  -  Xeq)c/d 


6Xeqc/d 


(2.6) 

(2.7) 


by  Eq.  (2.1). 

Eqs.  (2.1),  (2.2),  (2.5)  and  (2.7)  compose  a  phenomen¬ 
ologically  completedescription  of  the  transition  process  for 
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incorporation  into  the  flow  equations.  In  order  to  illustrate 
some  of  their  features,  note  first  that  d/c  in  Eq.  (2,7) 
plays  the  role  of  a  relaxation  time.  For  d  =  0.1  nan  and 
c  =  5  mm/n-sec,  this  becomes  t  =  a/c  ~  .02  psecs,  a  much 

25 

shorter  relaxation  time  than  was  reported  by  Novikov  et  al. 
and  one  which  would  play  little  role  in  shock  observations. 

The  equilibrium  curve,  Eq.  (2.5),  can  be  illustrated 
as  follows  for  iron:  Take 

N  =  (No/2)L1  +  tanh((P-PB)/AP  )j  (2.3) 

>*1  =  (vo/vl}  "  1 

P(v^)  =  1.667  +  3.4 

Av  =  v^-v2  =  *00596  cc/g 

vQ  =  .1275  cc/g 

Pm  =  .175  megabars 

2 

N  ,  or/ d  and  AP  to  be  varied 
o 

The  results  of  these  computations  are  shown  in  Fig.  26. 

The  parameters  for  each  curve  are  given  in  Table  2.1.  It  is 
apparent  that  the  equilibrium  curve  can  be  shifted  quite 
arbitrarily  in  the  transition  region  with  even  such  a  simple 
model  as  this.  Since  thermodynamic  calculations  give  almost 
horizontal  adiabat  and  Hugoniot  curves  in  the  mixed  phase  region, 
it  is  conceivable  that  careful  shock  measurements  in  the  mixed 
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Fig. 2  6 


P-V  Curves  Obtained  by  Varying  Parameters 
of  Martensitic  Model 
(See  Table  2.1) 
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Table  2.1 


Parameters 

for  Fig. 

2.1 

Curve 

AP 

No 

a/d2 

Number 

Mesabars 

1 

.015 

100 

.05 

2 

.01 

100 

.05 

3 

.02 

100 

.05 

4 

.015 

200 

.05 

5 

.015 

500 

.05 

6 

.015 

100 

.1 

7 

.015 

100 

.2 

phase  region  can  be  used  to  shed  light  on  deviations  from 
thermodynamic  equilibrium. 
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» o s ' n * c  r- An  equation  of  state  suitable  for  calculating  the  compression 
of  a  melting  solid  is  described.  Some  elementary  ideas  about  melting 
are  reviewed  and  some  standard  relations  between  P  and  T  in  the  melting 
region  are  described.  The  equation  of  state  and  melting  law  are  com¬ 
bined  in  a  program  for  calculating  the  Hugoniot  through  the  mixed  phase 
region.  Results  are  described  for  lead,  which  melts  at  a  shock  pressure 
of  about/ 400  kilobars  with  a  Kennedy  eauation  and  700  kilobars  for  a 
Siraorv-etfuation . 

The  Eyring  theory  for  equation  of  state  of  liquids  is  examined 
for  argon,  and  Hugoniot  curves  are  calculated.  Calculations  agree  with 
♦■he  most  dense  case  of  van  Thiel  and  Alder  to  13  kilobars,  then  depart 
ramatically  from  measured  values. 

—  The  theory  of  plastic  wave  propagation  in  two-dimensions  is 
discussed  and  calculations  of  allowed  directions  are  described.  These 
will  ultimately  be  of  use  in  discussing  the  reflection  of  obliquely 
incident  waves  in  an  elastic-plastic  medium. 

Some  of  the  basic  physical  mechanisms  in  solid-solid  phase  tran¬ 
sitions  are  reviewed  and  the  applicability  of  thermodynamics  to  such 
transitions  is  brought  into  question.  An  elementary  model  for  a  non¬ 
equilibrium  transition  in  iron  is  suggested  and  p-v  calculations  are  made 
for  several  values  of  the  parameters.  It  is  evident  that  no  conclusions 
about  the  time  dependence  of  the  a-e  transition  can  be  drawn  from  second 
state  shock  measurements,  although  it  may  be  possible  to  infer  useful 
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